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Abstract 

The finite element simulation of dynamic wetting phenomena, requiring the computation of flow in a 
domain confined by intersecting a liquid-fiuid free surface and a liquid-solid interface, with the three-phase 
contact line moving across the solid, is considered. For this class of fiows, different finite element method 
(FEM) implementations have been proposed in the literature and these are seen to produce apparently 
contradictory results. The purpose of this paper is to develop a robust framework for the FEM simulation of 
these flows and then, by performing numerical experiments, provide guidelines for future investigations. In 
the new framework, the boundary conditions on the liquid-solid interface are implemented in a methodolog- 
ically similar way to those on the free surface so that the equations at the contact line, where the interfaces 
meet, are applied without any ad-hoc alterations. The new implementation removes the need for complex 
rotations of the momentum equations, usually required to apply boundary conditions normal and tangent 
to the solid surface. The developed code allows the convergence of the solution to be studied as the spatial 
resolution of the computational mesh is varied over many orders of magnitude. This makes it possible to 
provide practical recommendations on the spatial resolution required by a numerical scheme for a given set 
of non-dimensional similarity parameters. Furthermore, one can examine various implementations used in 
the literature and evaluate their performance. Finally, it is shown how the framework may be generalized 
to account for additional physical effects, such as gradients in surface tensions. A user-friendly step-by-step 
guide specifying the entire implementation and allowing the reader to easily reproduce all presented results 
is provided in the Appendix. 
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1. Introduction 

Processes in which a liquid wets or dewets a solid surface are ubiquitous throughout industry and nature. 
They are known as dynamic wetting flows [DH]. Flows from this class occur, for example, when a liquid-gas 
meniscus propagates through a tube (Figure [T^) and during the spreading of liquid drops over solid surfaces 
(Figure [T]d) . In industrial applications, dynamic wetting flows are usually referred to as 'coating flows', 
emphasizing the technological objective to 'coat' the solid substrate with a liquid. One can distinguish 
between discrete wetting, as for ink-jet printing [3111], and continuous coating, where the liquid is deposited 
as a continuous film [ilin], e.g. using a hquid curtain falling onto a moving solid substrate (Figure [l];) . It is 
well established that controlling the dynamics of wetting is critical to the manufacture of a defect-free coating 
at an optimal speed [7l |8] . Manipulation of the wetting dynamics is often achieved by introducing additional 
physical effects which either alter the surface tension of an interface such as, for example, surfactant transport 
[21 [TO] or influence both the interfaces and the bulk [TT] . 

A defining feature of flows involving dynamic wetting is that, as confirmed by numerous experiments, 
e.g. [TOHl4j . the free surface meets the solid boundary at an angle 9, known as the 'contact angle', which is 
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Figure 1: Examples of dynamic wetting flows: a) propagation of a meniscus through a capillary tube, b) drop spreading over 
a solid surface, c) liquid curtain falling and coating a moving solid substrate. 

less than 180° (Figure[T]). As a result, mathematically, one faces a problem of finding a solution in a domain 
with a non-smooth boundary with different boundary conditions for the bulk equations of fluid mechanics 
on the two boundaries that form the contact angle [TH [TH] . 

From the modeling viewpoint, the difficulties are two-fold. Firstly, regardless of the contact angle value, 
there is no solution to the Navier-Stokes equations that would simultaneously satisfy all classical boundary 
conditions on the free surface and the solid boundary [HUT]. This nonexistence of what would be a 'classical 
solution' for dynamic wetting problems (in a weaker formulation also described in terms of a non-integrable 
shear-stress singularity) came to be known as the 'moving contact-line problem' [T51 [T5] . 

The second difficulty in the modelling of dynamic wetting phenomena is that the value of the contact 
angle, which provides a boundary condition for the equation describing the free-surface shape and hence is 
vital for finding the configuration of the flow domain, depends not only on the speed at which the liquid- 
fluid-solid contact line moves across the solid surface but, as experiments show, also on the flow field in the 
vicinity of this contact line |14lll9j . Therefore, in order to model dynamic wetting and be able to incorporate 
this effect, one has to accurately describe the flow near the moving contact line. 

Theoretical works that aim at modelling the process of dynamic wetting as a solution to a clear-cut closed 
mathematical problem (see [2] for a review) invariably preserve the contact angle, and hence a piecewise 
smooth boundary of the flow domain, and focus on formulating boundary conditions that would resolve the 
moving contact-line problem, i.e. allow for the existence of a solution (or, in a weaker formulation, make the 
shear stress integrable). In the vicinity of the moving contact line, the 'extra' physics that kicks in and has 
to be incorporated into the boundary conditions to resolve the problem makes it necessary to generalize the 
no-slip condition on the solid surface to allow for slip, i.e. a difference between the tangential to the solid 
component of the liquid's velocity and that of the solid. The simplest example is the Navier slip condition 
PU] in which slip is proportional to the tangential stress acting from the liquid on the liquid-solid interface. 

The complexity of free-boundary problems with moving contact lines makes it almost inevitable that 
numerical methods are required. If the free surface position and evolution are to be tracked explicitly and 
with high accuracy, it is advantageous to use an arbitrary Lagrangian-Eulerian (ALE) finite element method 
(FEM) in which computational nodes on the free surface evolve with it whilst those in the bulk move in some 
other specified manner that prevents the elements from degeneration. This method easily handles complex 
geometries, allows for a natural incorporation of boundary conditions formulated in terms of stresses and 
has been successfully used to describe a variety of free-boundary flows, e.g. PTH2^ . 

In applying numerical methods, including the FEM, to dynamic wetting flows, one encounters two main 
difficulties: (a) the characteristic length scale over which the classical boundary conditions have to be altered, 
i.e. in which the no-slip condition is violated, is much smaller than the scale associated with the bulk flow 
and (b) the boundary conditions that are applied on the two parts of the boundary that form the contact 
angle, i.e. on the free surface and the solid, are of different type. 

The disparity of length scales in dynamic wetting problems is usually so large that any attempt to 
resolve both scales using a numerical method with a fixed spatial resolution will inevitably produce a 
computationally intractable problem. Consequently, previous numerical implementations of dynamic wetting 
models split into two groups: those that explicitly attempted to resolve the smaller length scale using 
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specially designed meshes, e.g. ^51 [50^55] . and, the majority, where no such attempt was made. 

The consequence of each interface having different types of boundary conditions is that the contact 
line, where the two interfaces meet and the contact angle must be specified, does not fit naturally into the 
numerical scheme. In the FEM literature many different schemes have been proposed for the contact line 
implementation and, as discussed in detail in this work, many of them give different results for the same 
problem. 

A priori, it is unclear, both qualitatively and quantitatively, how the aforementioned difficulties affect 
the accuracy of a numerical scheme and how localized to the contact line any errors will be. To resolve these 
issues, we develop a computational framework which both fully captures the dynamics of slip and allows 
the contact line to be incorporated in a natural way. Then, this numerical tool is used to examine previous 
approaches and provide guidelines for future investigations. 

There is a large body of literature on the use of the FEM to simulate fluid flow ^34j , and it is well known 
how to accurately incorporate liquid-fluid free surfaces into this method, e.g. |21H23( I35j . In this paper, 
we complement these works by providing a detailed exposition, including practical recommendations, of a 
framework for incorporating dynamic wetting phenomena in a regular way into a FEM scheme. 

In ^ we give the equations of a basic, representative, model of dynamic wetting. The standard FEM 
approach to solving the governing equations is then outlined in Sj3j It is then described how previous 
works have attempted a number of different computational approaches to incorporate the dynamic wetting 
process. In fj4j we present a new framework which, by implementing the boundary conditions on the liquid- 
solid interface in a methodologically similar way to that used on the free surface, allows the contact line to 
be incorporated into the finite element procedure in a natural way. In the main body of the text the new 
approach is constructed so that it may be applied in an arbitrary coordinate system and is not dependent 
on the specific features of the FEM scheme, e.g. mesh design, which have been chosen. In the Appendix, 
we describe, in a step-by-step user-friendly guide, our FEM implementation in detail for the cases of two- 
dimensional and three-dimensional axisymmetric flow. 

In ^|5] we consider a representative example: the steady propagation of a liquid-gas meniscus through a 
capillary tube. This allows us to validate the new implementation against analytic results for spcciflc flow 
regimes and check the convergence of the code as the spatial resolution of the FEM mesh is increased. The 
code allows one to critically analyze previous numerical implementations of dynamic wetting models, com- 
pare them to the new one and provide practical quantitative recommendations for future code development. 

In ^ we use the new framework to show how 'additional' surface physics can be naturally incorporated 
into the FEM procedure and the effect which the dynamics has on the flow. Finally, in S|7j we summarize 
our findings and suggest future directions of research. 



2. A representative model of dynamic wetting flows 

Consider the steady spreading of a Newtonian liquid, with constant density p and viscosity ^, over a 
chemically homogeneous smooth solid surface. The liquid is surrounded by a gas which, for simplicity, is 
assumed to be inviscid and dynamically passive, of a constant pressure pg. Let the flow be characterized by 
scales for length L, velocity U, pressure fiU/L and external body force Fq. In the dimcnsionless form, the 
continuity and momentum balance equations are then given by 

V • u = 0, i?e u • Vu = V • P -I- S** F, (1) 

where 



-pi 



Vu+(Vu)^ , (2) 



is the stress tensor, u and p are the liquid's velocity and pressure, F is the external force density and I 
is the metric tensor of the coordinate system. The non-dimensional parameters are the Reynolds number 
Re = pUL/fjL and the Stokes number St = pFoL^ / (pU). 

Boundary conditions to the bulk equations are required at the liquid-gas free surface x = Xi(s, t), whose 
position is found as part of the solution, and at the liquid-solid interface x = X2(s,t), whose position is 
known; {s,t) are the coordinates that parameterize the surfaces. 
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On the free surface the kinematic condition of impermeability is 



u • ni = 0. (3) 

Hereafter, n is the unit normal to a surface pointing into the liquid, and subscripts 1 and 2 refer to the free 
surface and solid surface, respectively. 

In what follows, it is convenient to introduce a (symmetric) tensor I — nn, which is essentially a metric 
tensor on the surface. If an arbitrary vector a is decomposed into a scalar normal component a„ = a • n and 
a vector tangential part a| | (the subscript 1 1 henceforth denotes the tangential component of a vector) , so 
that a = a|| + a„n, we can see that, because n • (I — nn) = 0, the tensor (I — nn) extracts the component 
of a vector which is tangential to the surface, i.e. a • (I — nn) = a||. Then, = (I — nn) • V is the 
two-dimensional gradient operator playing the same role on the surface as V does in the bulk. 

The classical boundary conditions for the tangential and normal stress on the free surface can now be 
written down as 



ni 



Vu+(Vu)^ • (I - nini) + ^— ^ = 0, -p + ni ■ Vu+(Vu)' ■ ni = ' "\ (4) 

J Ca L J Ca 

Here, the superscript T denotes the transposition, cti is the surface tension of the interface and the capillary 
number Ca — ^U/aie is based on its equilibrium value cie. The pressure in the liquid is considered relative 
to its value Pg in the gas. 

The solid is considered rigid, impermeable and moving at velocity U, so that at the liquid-solid interface 
the standard form of the Navier-slip boundary condition, which in the present formulation replaces the 
condition of no-slip, is 

£^=^(u||-U||), (u-U).n2=0, (5) 

where /3 is the non-dimensional coefficient of sli]^ The no-slip condition is recovered in the limit f3~^ 0. 
The first term in ([5| represents the tangential stress acting on the interface and may be rewritten, to keep 
a consistent notation, in terms of P, so that conditions ^ take the form 

n2-P-(I-n2n2)=/3(u||-U||), (u - U) • n2 - 0. (6) 

The differential term ctiV • ni in the normal stress equation Q indicates that this equation requires its 
own boundary condition where the free surface terminates, i.e. at the contact line, where ni is prescribed 
provided that one specifies the contact angle 9 at which the free surface meets the solid surface: 

ni • n2 = —coaO. (7) 

Here, n2 is determined by the (known) shape of the solid, the contact angle 6 is specified by a particular 
model (or simply prescribed) and consequently ni may be determined and provides a boundary condition 
for the normal stress equation Q. 

Equations ([T]) with boundary conditions ([s]), Q, ([6| and ([t]) are the 'conventional' ones to describe flow 
near a moving contact line. Although the shortcomings of these conventional boundary conditions as a 
mathematical model for dynamic wetting are well known (see, for example [2]), as far as numerics is con- 
cerned, the described formulation is a necessary intermediate step towards implementing more sophisticated 
models. 

Next, we shall show how the FEM converts the outlined continuous problem into an approximate discrete 
one which may be solved computationally to a desired degree of accuracy. 



aiW ■ ni 



3. Conventional FEM simulation of dynamic wetting flows 

In this section, we will briefly recapitulate some of the main elements of the FEM, which is required in 
what follows as well as in the Appendix where a step-by-step implementation of the developed procedure is 
given. 



^Notc that in some works /9 is defined as the inverse of what we have chosen, e.g. |21l I36| 
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3.1. Weak form of the governing equations 

The defining feature of the FEM is that the computational domain V is tessellated into a finite number 
of non-overlapping elements, each containing a fixed number of nodes at which the functions' values are 
determined. Between these nodes the functions are approximated using interpolating functions whose func- 
tional dependence on position is chosen. Let Np be the total number of nodes in V at which the pressure 
is determined, be the number of nodes at which the velocity components are to be found and Ni be 
the number of nodes on the free surface Ai. The functions are approximated as a linear combination of 
interpolating functions each weighted by the corresponding nodal value. This gives a trial solution of the 
form 

Np Ar„ Ni 

p* = ^Pj?/'j(x), u* = ^Uj(/)j(x), X* = ^xij0ij(s,t), (xey, xi e Ai), (8) 
j=i j=i j=i 

where V'j, 0j and j are the interpolating functions and the asterisk is used to distinguish the approximate 
representations from the exact solution (p, u, Xi). The interpolating functions are constructed so that they 
take the value one at the node with which they are associated j and zero at all the other nodes to ensure 
that (pj,Uj,Xij) are nodal values. Note that here we use Roman letters for the indices (i,j, etc) to refer 
to the nodal values and approximating functions spanning the whole domain (globally); in the Appendix, 
where all the numerical details are given, these indices will be used alongside italicized ones (i, j, etc), 
which will refer to local, element-based, values and functions. It is not necessary to use the same level of 
approximation for the shape of the free surface and the velocity on it, but it is assumed here. In practice, 
each interpolating function is constructed in an element-by-element fashion and will only be non-zero over 
a handful of elements: this procedure is described in the Appendix but is not required for the following 
exposition. 

The representation (p*, u*, x^) of the actual solution (p, u, xi) is not exact and hence the substitution of 
([s]) into the governing equations ([T])-([2]) and the kinematic boundary condition Q will result in errors (the 
r's): 

V-u*=r'^, (9) 
e„ • [i?e u* • Vu* - V -P* - S't F] = r*^'", P* = -p*I + [vu* + (Vu*)^j , (10) 

u*-n^=r^ at x = x^, (11) 

where the unit basis vectors of the coordinate system Gq, {a — 1,2,3) have been used to obtain scalar 
equations. Hereafter, the superscript C will be used for quantities associated with the continuity of mass 
equation, M for the momentum equations and K for the kinematic boundary condition. We have chosen to 
use the kinematic condition on the free surface (|3| as the equation which determines the surface's position, 
as opposed to the normal stress condition Q which is shown in J 3.3 to be incorporated into the weak form 



of the momentum residuals. Then, given these three equations ([9|)-( 11 1 the problem becomes symmetric 
with respect to the boundary conditions: on each interface, i.e. the free surface and the solid boundary, two 
boundary conditions are applied. Henceforth, we shall drop the asterisks and realize that all functions to 
be considered represent approximate solutions to the problem. 

The aim is to minimize the errors spatially over the entire domain. In the method of weighted residuals 
[57] on which the FEM is based, this is achieved by minimizing a weighted spatial average of the error. 
Specifically, we multiply a given error by an arbitrary test function (the weight) W and integrate the equation 
over the corresponding domain. Applying this procedure to all errors creates the so-called residual equations. 
For the bulk equations this requires integration over V whilst for the kinematic equation integration is over 
the entire free surface Ai, so that the residuals to be considered are 



= [ W^{-K)r^ dV, R^^'" = f W'^^i^y^'" dV, R'^ = f W^{sA)r^ dA^. (12) 

Jv Jv J A, 

= = = 0, (13) 



IV JV JAi 

Insisting that these integrals vanish 
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for all test functions, of which there are infinitely many, creates the so-called weak form of the (strong form) 
governing equations of f|2j A link between the weak and strong forms is provided by the fundamental lemma 
of the calculus of variations which states that if certain equation-specific constraints on the test functions are 
satisfied, such as continuity and boundedness, and the relations above hold for all such test functions, then 
an approximate solution which is continuous will be exact, i.e. a solution of the strong form, see [38J p. 58] 
for specific details. In practice we choose a finite set of test functions to generate a discrete set of equations 
for the nodal values of the functions, which, it is expected, will approximate the continuous problem as 
closely as required. 

3.2. The Galerkin method 

It is useful to think of an equation associated with each node as 'determining' the nodal value of a 
function there. Mathematically, this is crude, but numerically it is helpful, and we henceforth use the word 
'determine' in this sense. Specifically, the Np continuity of mass equations can be thought of as determining 
the Np pressure nodal unknowns, the momentum equations in each coordinate direction determine the 
Nu velocity nodal unknowns in that direction and the A^i kinematic equations on the free surface determine 
the Ni nodal unknowns which specify the free surface shape (see the Appendix for a specific free surface 
representation) . 

To discretize the equations, one has to replace the infinite set of test functions, the W^s, by a finite set 
W ^ Wi {i — 1-N) which spans the solution space as completely as possible. In the Galerkin method 
this is achieved by choosing the test functions to be the same interpolating functions used for our variables. 
Specifically, the interpolating function of a given variable is the same as the test function for the equation 
which determines that variable's value, so that 

Wi'^(x) = Vi(x) (i=l-iVp), Wi*^(x) = (/)i(x) (i-l-A^„), Wi^(s,t) = (/.i.i(s,t) (i = l-A^i). 

(14) 



After substituting the interpolating functions in (14 1 into the expressions for residuals (12), equations n3\ 



are transformed into a finite set of algebraic equations 

i?P = (i=l-7Vp), R^^'^^O {i^l~Nu), R^ = Q (i=l-iVi). (15) 

In |3.3[ we give specific expressions for these and show how they can be manipulated to allow some of the 
boundary conditions to be easily incorporated. As one increases the number of nodes in the domain, the 
error from approximating a continuous problem with a discrete one should decrease, and this will be checked 
for the problem in question in §5.3| 

Here, global residuals and functions, i.e. quantities considered over the entire domain, have been con- 
sidered. In practice, the global interpolating functions are constructed piecewise over elements, i.e. locally, 
as low order polynomials. Then, the global residuals are obtained by summing up contributions from each 
element. In the velocity-pressure formulation of the incompressible Navier-Stokes equations, mixed interpo- 
lation [35], with velocity approximated by a polynomial of degree at least one more than pressure, ensures 
that the Ladyzhenskaya-Babuska-Brezzi [lOj condition is satisfiecQ Instead of considering the elements of 
varying shapes in the global domain, it is useful to map each element onto a master element over which 
analysis can be systematically performed: this ensures that the solution to problems in complex domains 
pose no additional difficulty. This procedure is described in detail in the Appendix, but will not be required 
in what follows. 

3.3. Discrete problem formulation 



In this section, using ( 12 ) and ( 14 1 the discretized weak form of the governing equations ( 15 ) is pre- 
sented and it is shown how the boundary conditions are applied. Functions appearing in the residuals are 
approximated using (|8|. 



* Equal-order methods may also be used with stabilization, e.g. I41II42I . but these issues lie beyond the scope of this paper. 
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The continuity of mass (incompressibility of the fluid) residuals i?p are 

r9 = I ^Pi\/ -udV (i = l~Np), (16) 
Jv 



Momentum equations are 



/ <?!)ie„ • [i?e u • Vu - V • P - S** F] dV" (i = 1~7V„). (17) 
Jv 



Assuming that the interpolating function is at least once differentiable, the stress tensor's contribution can 
be represented as 

^ 0ie„ • V • P = -V • (0ie„ • P) + V (0ie„) : P, (18) 
so that, after applying the divergence theorem, one has 

R^^'" = I [(ji-^e^ • (i?e u • Vu - St F) + V(0ie„) : P] dV + [ (ji-.e^ • P • n dA, (19) 

Jv J A 

where, as before, n is the inward normal to the surface A which forms the boundary to V . Therefore, the 
volume contribution to the momentum residuals is 

{R-^ "^ ^ = j [(/)ie„ • (i?e u • Vu - 5i • F) + V(0ie„) : P] dV, (20) 
whilst the contribution to the momentum residual from the surface is 



!)ie„-P-ndA. (21) 



In the last expression, only when node i is on the surface A will be non-zero, i.e. only nodes on the surface 



of V will give contributions to the momentum residuals via (21). 

Previously, it was described how for each nodal unknown there is an equation which is assumed to de- 
termine its value. In all numerical methods, to maintain the equality between the number of unknowns 
and equations, in applying the boundary conditions we must either (a) replace a bulk equation at boundary 
nodes by an equation representing the discretized form of the boundary condition, or (b) manipulate the 
bulk equations to allow the boundary condition to be incorporated into them. We refer to (a) as an essential 
imposition of the boundary condition; the bulk equation which is replaced is referred to as the dropped equa- 
tion. The (b)-type implementation of boundary conditions will be referred to as natural imposition, which. 



for stress conditions, is made possible by the contribution (21 1 that is replaced by the corresponding expres- 
sion for the stress on a given surface (the FEM is therefore referred to as allowing a natural incorporation 
of stress-type boundary conditions). 

On the free surface, boundary conditions ([4]) can be combined into a more computationally favourable 
form 

Ca ni • P + V" • [cTi (I - nini)] = 0, (22) 
where cti (I — riirii) is the surface stress, playing the same role on the surface as P does in the bulk, so 



that ( 22 ) is actually a balance of stress acting on the interface (first term) and that acting in the interface 



(second term). Then, using (22), one can rewrite (21) on the free surface as 

/ 0ie„ -P-m dA = --^ / 0i.ie„-V"-[cri(I-nini)] dAi. (23) 
Jai (^o,Jai 



Direct calculation of the integrand in ( 23 ) would require the evaluation of the derivatives of the vector normal 
to the surface. It was initially suggested by Ruschak |43], and generalized for three-dimensional problems 
in |23[ I44j , that by using the surface divergence theorem one could obtain an expression in which no such 
derivatives are required. Lowering the highest derivatives reduces the constraints on the differentiability of 
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the interpolating functions which are used to approximate the free surface shape, i.e. (^i^i, and hence aUows 
lower-order polynomials to be used. The required transformation of ( 23 ) is given by 



/ 0ie„ • P • ni dAi = / [V* • (0i,iCTie„ • (I - riirii)) - ctiV" • {(f>i.iea)] dAi. 



(24) 



Then, using the surface divergence theorem [151 p. 224], which for a surface vector a with no normal 
component a = a| | is 

V"-a|| dA^- a|| mdC, (25) 

A Jc 

where the unit vector m is the inwardly facing normal to the contour C that confines A (Figure [2]), with 
a|| = (pi^iaiBa ■ (I — nirii) we obtain 

iBa • P • ni dAi = ^ I cTiV • {(l)i.\ea) dAi + [ 0i,icrie„ • mi dCi. (26) 
Ai L^aJAi (^a Jc^ 



Thus, on the free surface, the term (21) in the momentum residual is now split into surface and line 
contributions 



= 7^ / CTiV* • (^i^iGa) dAi, f-R^'") = 7^ / 0i,iCTiea • mi dCi, 
\ ' J Ai Ca Ja, V Jci CaJc, 



(27) 



with the stress boundary conditions (22) incorporated. 

Notably, the same procedure of integrating by parts and using the divergence theorem has been used on 
both the surface stress term • [ai (I — nn)] and the bulk stress term V • P. In both cases, this has created 
contributions from the boundary of that term's domain, i.e. the confining contour and surface, respectively. 
Consequently, the momentum residual now contains a cascade of scales 

M,a _ ( r>M,a\ . ( jyM ,a\ . ( r>M ,a 



K'" = i< (28) 

which represent, respectively, the volume, surface and line contributions. In particular, part of the contour 
which bounds the free surface is the contact line Cd where the liquid-gas interface meets the solid. Other 
boundaries to the free surface further away from the contact line, for example an axis of symmetry, are not 
considered here. 

At the contact line, it is useful to rearrange the term in the integrand of the contour integral in (27) 
by representing the vector mi in terms of a linear combination of its components parallel to n2 and m.2 
(Figure [2]): 

mi = (mi • m2) m2 + (mi • n2) n2. (29) 

Using this to impose the contact angle ([T]), which is a boundary condition for the free surface shape, results 
in 

[Rf'") = ai4>iea ■ {ni2Cos9 + n2sm6) dCd, (30) 

V /Co! oa Jq^^ 

where the normal and tangent vectors on the solid surface are known a-priori. If contact line's tangent 
vector is t^, then m2 — it^ x n2, with the sign chosen to ensure the inward facing vector m2 is selected. 

Thus, the contact angle equation ([t]) can be applied in a natural way, that is without needing to drop 
another equation in order to impose the contact angle condition and thus fix the shape of the free surface 
at the contact line. 

The kinematic condition residuals Rf^ complete the residuals on the free surface: 

Rf ^ I 0i,iU • ni d^i (i=l-7Vi). (31) 

JAi 

Thus far, the implementation has been straightforward: the stress conditions are applied in the weak 
form of the momentum equation, and their integration by parts allows specification of the angle at which 
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Figure 2: Illustration showing the vectors associated with the liquid-gas free surface Ai and the liquid-solid surface A2 in the 
vicinity of the contact line Cd- 



the free surface meets the sohd boundary. This means that momentum equations on the free surface do not 
need to be replaced in order to apply essential surface or line conditions. In other words, all of the bulk 
equations are applied at every node, to determine the velocity components, with both the free surface stress 
conditions and the specification of the contact angle fitting naturally into the weak form of the momentum 
equations. The kinematic condition is applied as a separate equation to determine the free surface position. 
Notably, the implementation is independent of the free surface shape, that is the curved nature of the free 
surface is as easy to handle as, say, a planar surface which aligns with coordinate axes. 

On the rigid solid surface, where we have A'2 nodes, it is the Navier condition and the condition of 
impermeability (|6| that have to be imposed. The Navier-slip condition is of stress type and hence contributes 
to the weak form of the momentum equations in a natural way, giving 

'")^^ = 02,iea • (u|| - U||) dA2, (32) 

where 02.1 is an interpolating function for the liquid-solid interface corresponding to the i-th node. However, 
the impermeability condition is not of stress type and hence does not fit naturally into the weak form of 
the momentum equations. Consequently, it is usually applied as an essential condition which creates a new 
residual equation at each of the nodes on the liquid-solid interface: 

R( = (i = I-N2) where R( = / 02.1 (u - U) • n2 dA2. (33) 

Ja2 

Then, as previously described, in order to retain equality between the number of unknowns and number 



of equations in the numerical scheme, an equation must be dropped to make room for ( 33 ) . As ( 33 1 is an 
equation which determines a component of velocity, and it is the momentum equations that are regarded 
as determining these, it is a momentum equation that should be dropped. The dropping of a momentum 
equation in order to make space for an essential condition is achieved more formally, see e.g. |34[ 138] . by 
defining the interpolating functions such that they are zero at nodes at which an essential condition is to be 
applied. Either way, the end result is exactly the same. 

The procedure outlined above is straightforwardly achieved when the solid surface is parallel to a coor- 



dinate plane, so that one can drop the momentum equation normal to this plane and apply (33 1 instead. 
The momentum equations tangential to the boundary only require an expression for the tangential stress 
and this is provided by the Navier-slip condition. If the solid surface is not parallel to a coordinate plane, 
then, so as not to favour one of the axes, before applying the boundary conditions one must project the 
momentum equations at each node on the surface so that they are tangent and normal to it |46j . This is 
not a complex procedure, but it is considerably more cumbersome, and hence more time-consuming, than 
that encountered on the free surface, where curved shapes do not cause any additional difficulties. One can 
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imagine that such a procedure will be computationally intricate when considering flows over topographically 
complex surfaces where the momentum residuals will have to be projected on different axes at each node. 

So far, the weak formulation of the bulk equations has been derived and it has been shown how the 
boundary conditions on the liquid-gas free surface and the liquid-solid interface are implemented. Addi- 
tionally, it was shown that the contact angle may be naturally imposed by specifying a contribution to the 
momentum equations from the contact line. So, on the liquid-gas free surface all momentum equations are 
applied, with stress-type boundary conditions naturally imposed, as well as the kinematic condition, whereas 
on the liquid-sold interface one momentum equation was dropped in order to apply the essential condition of 
impermeability. Next, we consider what happens when we reach the nodes at which these interfaces, which, 
as described above, are treated in methodologically different ways, meet at the contact line. 



3.4- Conventional FEM implementation of the contact line 

Ideally, one would use an implementation which allows all the boundary conditions from both sides of 
the contact line to be applied in addition to the contact angle specification at the contact line. However, 
it will be shown that a count of the number of equations and unknowns indicates that this is not possible. 
Consequently, in the method as it stands, one has to favour certain boundary conditions over others to 
arrive at an equality between equations and unknowns. 

To make the argument most transparent, consider as an example the free surface meeting a solid which is 
parallel to the z-axis and motion occurs in the (r, z)-plane of a Cartesian coordinate system with the liquid 
below, in terms of z, the free surface (e.g. see the geometry of Figure |4|. In this geometry, the free surface 
can be represented, for example, as z = h{s), where s is the arclength along it, so that in the discretized 
form hi (i = l-A'^i) are the unknowns that determine the free surface's position. 

The incompressibility equation is applied at every node and determines the pressure. At every node of 
the liquid-gas free surface both momentum equations are used, with the stress boundary conditions applied 
naturally, for the components of velocity (ui, Wi), and one kinematic condition is applied for the free surface 
unknown hi. On the liquid-solid interface the ^-momentum equation is applied to determine Wi, with the 
Navier boundary condition applied naturally, whilst the r-momentum equation is dropped to allow the 
essential impermeability condition, Ui = 0, to be satisfied. 

At the contact line, node i = 1, we have unknowns ui,wi,hi, with hi corresponding to the z-coordinate 
of the contact line. Application of the z- momentum equation, to determine wi, allows both the Navier 
boundary condition along the solid and the z-projection of the stress conditions on the free surface to be 
imposed. This is a good start, but there are now three equations left - the r-momentum equation, the 
kinematic condition on the free surface and the impermeability condition on the solid surface - for only two 
unknowns ui,hi. The conventional solution to this difficulty is to favour the latter two equations and not 
apply the r-momentum equation at the contact line. This is in keeping with the usual FEM approach of 
replacing bulk equations wherever essential boundary are applied. However, the equations at the contact 
line now enforce the z-projection of the free surface stress boundary condition but not the r-projection. 
Furthermore, the contact angle specification is no longer applied into the r-momentum equation at the 



contact line. Using (30) and noting that m2 = (0,-1) and n2 = (—1,0), the contribution from the contour 



integral in ( 30 ) to the z-momentum equation at the contact line will be 



(ijf'^) = -01,1(71 cos 6'/Ca. (34) 

In §3.5[ implementations of the contact line in FEM codes known in the literature are discussed. In 
all these cases, as above, at least one momentum equation is not applied at the contact line. This type 
of implementation, in which equations are dropped at the contact line, will be referred to as Type A and 
in this framework the following two different implementations of the contact angle specification have been 
developed: 



(A-a) The contact angle is imposed naturally, as described above, using ([30| (whilst one momentum equation 
is dropped already). 
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Figure 3: Illustration of a solution in which the contact angle at which the computed free surface meets the solid dc differs 
from the applied angle B which appears in the mathematical model. 

(A-b) The contact angle is imposed as an essential condition (with one more momentum equation also 
dropped). 

Exactly the same methodological issues arise in fully three-dimensional flow simulations. 

Before reviewing previous works that use the FEM to simulate dynamic wetting flows, with particular 
attention to the development of different implementations of the contact angle, it is necessary to make a 
remark concerning the different 'contact angles' now involved in the modelling. When the contact angle 
specified in the mathematical model is incorporated into the FEM code in a natural way, i.e. absorbed 
into the contribution to the residual from the contact line, as in the (A-a) implementation, there is no 
guarantee that the contact angle that emerges after computations are performed, i.e. the angle at which the 
computed free surface meets the solid boundary, will be equal to the specified one. We will refer to the angle 
incorporated naturally into the code, and appearing in the mathematical model, as applied, and denote it 9, 
whereas the contact angle at which the computed free surface meets the solid will be referred to as computed 
and denoted as 9^ (Figure [s]). In contrast, in the (A-b) implementation there is a single equation which 
ensures that 9c — 9. 

3.5. Previous approaches 

We will now consider the development of the two different FEM implementations, (A-a) and (A-b) in 
the literature and, in particular, focus on how they appear to produce contradictory results. This is by no 
means supposed to be an exhaustive review, rather a pointer to some key papers and their conclusions. 

Building on early articles considering the FEM simulation of free surface flows, for example |37], pio- 
neering research into the use of FEM to describe the dynamics of coating flows, and hence dynamic wetting 
phenomena, came from the University of Minnesota. The study of flows with contact lines began by com- 
puting static meniscus shapes in I48-50| after which problems such as the die-swell phenomenon, with static 
contact lines, were simulated [Sir 53 . and coating processes involving dynamic contact lines considered in 
pn 1221 [511 [5S]. The conclusions of these investigations and the techniques, subsequently used in many 
works, are as follows. 

Initially, the Minnesota group used the weak formulation to impose the shape of a free surface where 
it intersected a boundary, i.e. implementation (A-a). This appears to have worked well for imposing the 
contact angle of a static meniscus shape. Using this implementation to impose the contact angle in a dynamic 
wetting problem was discussed in ^54; . p. 263]. However, in |21l p. 217], it was found that implementation 
(A-a) used in [53] resulted in the computed angle 9c differing significantly from the applied angle 9, and this 
discrepancy, in the worst case, even prevented the numerical scheme from obtaining a solution. To circumvent 
the problems observed in j21l . it was suggested that all the momentum equations should be dropped at the 
contact line to allow the contact angle to be imposed as an essential condition, i.e. implementation (A-b) was 
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proposed. Similar conclusions are discussed in 55', p. 333] and |221 p. 46], and it appears that this method 
was henceforth adopted. These papers appear to be the first in which it is recognized that implementation 
(A-a) does not guarantee that the applied angle will be equal to the computed one. Their conclusion is that 
implementation (A-a) must be utilized to guarantee that 9^ = and allow seemingly satisfactory results to 
be obtained. 

In [56J , an FEM code is designed for the simulation of free surface flows and implementation (A-a) is used 
to specify the contact angle. This paper is the first to suggest a new interpretation of the phenomenon which 
the Minnesota group had observed. The authors stress that (A-a) does not, and should not necessarily, make 
Oc equal to 6: the term (30), where the contact angle is incorporated into the weak formulation, merely 
creates a net interfacial force on the contact line. The computed angle O^, which varies from the applied one 
0, is then referred to as the 'dynamic' contact angle, which, for small capillary numbers, was found to be close 
to 9. This interpretation was adopted in [571 158j . where implementation (A-a) was used in studies of drop 
impact and spreading phenomena. As the drop spreading process commences, with the initial computed 
contact angle at 9c = 180°, for early times there is a huge, unavoidable, difference between 9 (in the model 
it is assumed that is a constant which in all cases is less than 100°) and 9c- From these publications the 
idea emerged that the dynamic behaviour of the computed angle, which differs from the (constant) applied 
one, is a physical effect, despite the fact that there is no 9c in the original mathematical formulation of the 
problem whose solution the code is supposed to describe. 

In all of the aforementioned papers, the FEM mesh has been designed to capture the general character- 
istics of the bulk flow, and the region near the contact line is given no special treatment. In particular, there 
is no attempt to capture the dynamics of slip, which approximately occurs on the (non-dimensional) slip 
length , by placing a number of elements within this distance of the contact line sufficient to resolve this 
region. The articles use what is referred to as 'numerical slip', which broadly refers to numerical schemes 
where slip is only prescribed, often in a couple of elements adjacent to the contact line, to ensure that the 
nodes at the contact line are not pinned to the solid by the no-slip condition. In other words, slip is used 
as a numerical convenience and no attention is paid to the extra physics which it introduces that occurs 
on a length scale greatly smaller than the elements used. Given that the region of numerical slip depends 
on the size of the elements adjacent to the contact line, the results of such computations are obviously 
mesh-dependent. 

An alternative approach is to create a special mesh which allows small elements to be placed around the 
contact line and the slip dynamics to be captured. This is achieved in [351 [33] where the propagation of a 
liquid-fluid meniscus through a capillary tube is considered. In 133!, implementation (A-a) is used, whilst in 
[32j the normal-stress boundary condition is used to iterate the free surface position with the contact angle 
applied as an integral constraint on the pressure level, i.e. it is not imposed as an essential condition. In [33], 
it is reported that the computed angle 9c was within 0.01° of 9 and in |32j no significant deviation is noted. 
These results appear to suggest that, if sufficiently small elements are used near the contact line, so that 
the dynamics of slip is well resolved, then the weak formulation can indeed be used to impose a specified 
contact angle. 

Even in recent papers, there is still much disagreement over which approach is correct. For example, 
in |59j approach (A-a) is used, with the computed angle, which differs from the applied one, interpreted as 
the dynamic contact angle; whereas, say, in [36] the implementation (A-b) is taken, with all the momentum 
equations dropped at the contact line. Variance also exists in whether or not the dynamics of slip is resolved: 
some authors design special meshes to do so, e.g. in [23131] [SD], whilst in many other publications the region 
of slip is under-resolved either because of constraints on computational power, as could be the case for 
three-dimensional simulations, or simply because this phenomenon is not considered important. 

Previous investigations show that when a 'standard' mesh is used, i.e. one which does not incorporate 
smaller elements near the contact line, implementations (A-a) and (A-b) certainly create different solutions: 
in case (A-a) it has been established that 9c deviates from 9 imposed via the weak formulation, whereas 
in (A-b) the contact angle is imposed as an essential condition at the expense of dropping another of the 
momentum equations at the contact line, so that it is unclear whether or not the imposition of 9c — 9 creates 
problems with the accuracy elsewhere in the scheme. The deviation of 9c from 9 which was observed when 
using (A-a) has been attributed to the fact that the implementation only imposes a net interfacial force at 
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the contact line, but, as far as we are aware, there has been no thorough convergence studies performed to 
see if the contact angle varies as the mesh is refined and whether additional mesh fineness in the immediate 
vicinity of the contact line has an influence on the flow further away. 

What has been presented thus far broadly summarizes the popular implementations used in the literature. 
However, a closer look at previous works reveals even more sub-varieties of implementation. For example, 
in implementation (A-a) it is usually the projection of the momentum equation tangent to the solid surface 
which is retained; however, in |54j the projection of the momentum equation onto the free surface's normal 
is chosen. Furthermore, in j60j . instead of applying the impermeability equation separately on the free 
surface and the solid surface at the contact line, i.e. creating two FEM residual equations, the two were 
combined to form one impermeability residual with contributions from the two surfaces. Encouragingly, this 
created room in the formulation to apply both momentum equations at the contact line alongside the (one) 
impermeability equation. However, it was observed that using this approach resulted in an 0(1) velocity 
at the contact line, which violated the requirement that both the free surface and the solid surface were 
impermeable. Implementation (A-b) was then adopted to overcome this issue. No doubt, other variations 
exist and clearly the implementation options at the contact line are plentiful. 

By now, it is unclear as to which implementation of the above, if any, is correct, what effect the dominant 
parameters in the contact line implementation (3 and Ca have on the seemingly contradictory outcomes and 
how the spatial resolution of the mesh influences the results. In the next section, we describe an approach 
which ensures that both the liquid-solid and liquid-gas interfaces are treated in a methodologically similar 
way so that the contact line does not present an obstacle to the numerical implementation. This approach 
does not contain any loose-ends at the contact line and hence does not allow for creativity on the part of 
the user. After rigorously validating our numerical code we will use it to address the aforementioned issues 
and provide practical advice for future code development. 

4. New approach 

In the two implementations (A-a) and (A-b) that have been described and which are both frequently used 
in the FEM literature, at least one of the momentum equations is dropped at the contact line. Consequently, 
some of the naturally imposed stress-type boundary conditions and line contributions are dropped. By now, 
it is unclear whether these droppings at the contact line will cause the numerical scheme to produce a 
spurious solution or whether there will be an influence on the speed of convergence of the code towards a 
correct solution. To see if this is the case, we develop a new implementation which requires no equations to 
be dropped at the contact line and thus allows us to examine how dropping equations at the contact line 
influences accuracy. To achieve this, we will outline a known, but rarely used, method for incorporating 
boundary conditions on the liquid-solid interface and show, for the first time, how this gives a perfect 
framework for incorporating the contact line. 

It was shown in §3.3| that the stress boundary conditions on the free surface are incorporated naturally 
into the weak form of the momentum equations. The kinematic condition ^ is then the extra equation to 
determine the free surface shape. However, because the solid surface's shape is prescribed, on the liquid-solid 
interface the normal stress is determined as part of the solution, as opposed to appearing in a boundary 
condition which would slot into the weak formulation. Consequently, there appears one less unknown on the 
liquid-solid interface than on the liquid-gas interface and hence an equation needs to be dropped to enable 
the essential boundary condition of impermeability of the solid surface equation to be satisfied. In other 
words, when the surface's shape is prescribed, the normal stress condition can no longer be applied and 
the kinematic condition plays the role of the second boundary condition for the Navier-Stokes equations. 
Clearly then, because the free surface and the solid surface are handled in different ways, the contact line 
where they join will have to be specially treated. 

In order to eliminate the asymmetry between the two interfaces, an extra unknown, denoted A, is 
introduced into the scheme; this unknown will be equal to the normal stress at the liquid-solid interface 

A = n2-P-n2. (35) 
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The introduction of A as an extra unknown makes it possible to reinstate the momentum equation normal 
to the solid surface at all nodes on the liquid-solid interface: all momentum equations will then be applied 
everywhere. In other words, on the free surface all momentum equations are applied with the impermeability 
condition used to determine the 'extra' unknowns specifying the position of the free surface, and now we 
have that all the momentum equations are also applied on the solid surface where the geometric constraint 
of the prescribed shape allows the impermeability condition to be used to determine the 'extra' unknown A, 
i.e. the normal stress. As a result, the free surface and the solid boundary appear to have been treated in 
the same way and hence, once the two boundaries meet at the contact line, the contact line conditions can 
be implemented naturally, without dropping any of the equations there. 

The new unknown is non-dimensionalized using the characteristic scale ^U/L and then approximated 
using the velocity interpolating function on the liquid-solid interface, so that 

N2 

A = £Aj02j. (36) 

Then, after splitting the stress into components normal and tangential to the solid 

e„ • P = e„ • [(n2 • P • na) na + na • P • (I - nana)] , (37) 
and applying the Navier boundary condition (|6| as before, the contribution from the liquid-solid interface 



to the momentum residuals (21 ) becomes 



02,i[A (e„.n2)+^e„.(u||-U||)] dA2, (38) 

with the first and second terms in the integrand associated with the normal and tangential stresses, respec- 
tively. Now, on both interfaces all momentum equations are applied to determine the velocity components 
and the extra kinematic-type equation determines the additional unknowns xi or A. 

The new approach introduces a significant simplification, at the expense of relatively few additional 
unknowns, when the solid surface is not parallel to a coordinate axis: instead of having to project the 
momentum equations in the normal and tangential directions to the surface, to allow the impermeability 
condition to be applied, no additional work is required as all equations are applied everywhere in a coordinate- 
free manner. This is particularly advantageous for the solid surface with complex topography. Furthermore, 
in contrast to implementations of Type A, all momentum equations are applied at the contact line, with 
contributions to the momentum residuals from each interface, so that the contact line equations no longer 
require special attention. In other words, the new approach allows the contact line to naturally fuse the 
liquid-solid and liquid-gas interfacial equations so that there is no scope for options on how the equations are 
implemented at the contact line, i.e. there are no 'loose ends'. This will be referred to as an implementation 
of Type B: 

(B) All momentum equations are applied at the contact line and the contact angle is imposed naturally 



into the weak form using ( 30 1 



The realization that this particular implementation of the boundary conditions on the liquid-solid in- 
terface allows the contact line to be handled naturally is new; however, the idea of introducing additional 
unknowns on the liquid-solid interface has been well documented. Building on the ideas of Babuska j61j . 
the method was considered by Verfiirth [52] for implementing slip conditions, and it was shown that A is the 
Lagrange multiplier of the impermeability equation (|6|. This is analogous to the situation with pressure, 
which is well known to be the Lagrange multiplier of the incompressibility constraint, as, in both cases, 
modelling assumptions (zero Mach number/rigid solid) result in a variable (p/A) being determined by an 
equation in which it does not appear (incompressibility/impermeability). In I63j, it is shown, for a case in 
which there are no contact lines, that results obtained from the approach of dropping the normal momen- 
tum and replacing it with the impermeability condition converge to the same computational results as the 
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Lagrange multiplier approach. Lagrange multipliers are also used in [64 to impose both the liquid-liquid 
free surface and the far-field boundary conditions. In [65], for the problem of flow over a backward facing 
step on which a slip boundary condition is applied, it is shown that the best results are produced when the 
Lagrange multiplier is approximated with the same basis functions as the velocity, as we have chosen to. 

The method described above ensures that contact lines may be treated in a methodologically consistent 
way with no equations dropped and hence no loose ends that leave room for user intervention and can 



lead to the myriad of approaches described in ^3.5 Now we shall investigate the accuracy of the new 



implementation, and then compare its results to those obtained using previous implementations (A-a) and 
(A-b). 



5. Numerical experiments: meniscus in a capillary tube 

The new implementation, as well as the ones considered in previous works, have been incorporated into 
the FEM scheme described in the Appendix. First, results from the new implementation are compared to 
those obtained from asymptotic analysis. Having confirmed the accuracy of the scheme, a practical guide 
to the spatial resolution which numerical schemes modelling dynamic wetting phenomena require is given 
in §5.3| after which the previous approaches are examined in §5.4| 

As a representative dynamic wetting flow, consider the steady propagation of a liquid-gas meniscus 
through a cylindrical capillary tube. This flow configuration has the advantage that the geometry is simple 
and that published results, which agree with each other, exist from spatially well resolved FEM schemes 

El Eg. 

Axisymmetric cylindrical polar coordinates (r, z) are used with the z-axis along the centre of the tube and 
pointing towards the gas; azimuthal symmetry about this axis is then assumed (Figure]!]). The characteristic 
length L and velocity scale U are provided by the radius of the capillary and the speed at which the meniscus 
propagates, respectively. The solution is computed in a frame in which the contact line is stationary, so that 
in this frame the solid, which is located at r = 1, moves with velocity U = — e^. 

The bulk equations arc (]!]), with the standard free-surface conditions (]3])-(]4]), the Navier-slip and im- 



permeability conditions (6 1 and new normal stress specification (35) at the liquid-solid interface, and the 
contact angle condition (7), with a prescribed velocity- independent value of 0, at the contact line whose 
location is (1,0). The surface tension of the free surface takes its constant equilibrium value of cti = 1. 
Additional boundary conditions at the axis of symmetry r — are 

u • e,. e,. • P • 6;, 0, (39) 

and the condition of smoothness on the free surface shape at r = that ni • e^ = 0. Boundary conditions 
in the truncated far field at z — —Zf are provided from the assumption that the flow is fully developed 

u^O, w^ar^ + b; a = -2/(4/3-^ + 1), 6 = -a/2; (40) 

where (u, w) are the components of velocity in the axisymmetric coordinate system, and if we have no-slip 
/3~^ = the usual Poiseuille profile is recovered. For a similar problem, in 03] it is noted that extending the 
truncated far-field further than two and a half radii away from the contact line did not alter their results. 
As the main interest here lies in the computational aspects near the contact line anyway, the truncated 
far-field is placed at zt— 3 and the issue is not considered further. The implementation of the boundary 



conditions ( 39 ) and ( 40 ) into the finite element procedure is carried out in the conventional way with the 
stress conditions incorporated naturally into the weak formulation and the Dirichlet conditions applied 
instead of momentum equations as essential conditions. New variables could have been introduced on these 
boundaries to ensure the momentum equations are always applied, but the implementation of the conditions 
on these boundaries into the FEM is not the focus of this paper. 

5.1. Parameter regime of interest 

To obtain estimates on the values of the similarity parameters, consider the displacement of air by water 
with viscosity /i = 1 mPa s and density p = 10'^ kg m~'^ in a capillary of radius 100 /im. The water-air 
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Figure 4: The computational domain for a liquid-gas meniscus propagating through a capillary tube with a close-up of the 
contact line region showing the coordinate system used in the derivation of the local asymptotic results. 

meniscus advances at a speed U = 1 m s^^ and has a surface tension of aie = 70 mN m^^. Finally, an 
estimate for /3 often used in the literature for Newtonian fluids is that /? ^ 10^/i. This gives 

that Re w 100, Ca w 10^^, f3 w 10^ and external forces are neglected so that St — 0. For simplicity and 
transparency, we consider acute contact angles to avoid the spurious numerical behaviour encountered, and 
then resolved, in ^67 for obtuse angles. In particular, 9 = 30° is taken as a base state. 

Initially, we shall consider Stokes flow, i.e. Re — 0, to allow for a comparison with asymptotic results. 
After confirming the accuracy of our approach, we consider more numerically challenging flows where the 
Reynolds number is non-zero and the capillary number is such that the free-surface becomes heavily de- 
formed. 



5.2. Comparison to local asymptotics 

To confirm the accuracy of our numerical results it will be useful to compare them to local asymptotics 
derived in the limit p — >■ where (p, 0) are polar coordinates centred at the contact line with the angular 
coordinate measured from the solid surface through the liquid (Figure |4| . Such an analysis can be found, 
for example, in [17] . The leading-order term of the streamfunction, introduced using 

1 (9V 

is given by 

7/, = p2^(e), F(e) = (Bi +536 + 53 sin 26 + ^4 cos 29), (42) 

where the constants of integration Bi {i = 1-4) are: Bi = —B4 = — /3/4, B2 = —Bi/9, B^ = Bi cot 26*. For 
comparison with our numerical results, we will use the radial velocity Up and the normal stress A on the 
solid surface = 0, which are given by 

Up^pF'ie^O); A = -4B2lnp-2(B2 + 2S3)+po, (43) 

where po is a constant of integration. Later, the free-surface shape will also be needed, which is parameterized 
(Figure [4]) as Q — Q{p), and locally has the form 

e ^ e + Ca/3[Ap\iip + Bp]+o{p), (44) 
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Figure 5: Left: streamlines in increments of i/i = —0.02 computed for Ca = 0.01, /3 = 10"", S = 30°. Right: corresponding mesh 
used in this region. 
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Figure 6: Left: streamlines near the contact line in increments of i/i = —2 X 10 ^ computed for Ca = 0.01, /3 = 10^,9 = 30°. 
Right: corresponding mesh used in this region showing a high degree of spatial resolution. 



where A^B are 0(1) functions of 0. Expressions in (43) will be sufficient to validate the numerical results. 
Note that to leading order the domain is wedge-shaped so that the radial velocity Up should closely ap- 
proximate the velocity tangential to the free surface Ug from the computed solution, which will be used for 
comparison. 

Consider typical parameter values Ca = 0.01, /3 = 10^, 9 — 30°. Figure [s] shows the corresponding 
streamlines and the mesh which was used for the computation. It can be seen that the motion of the solid 
causes fluid nearby to be dragged down and consequently, to ensure continuity of mass, this is replenished 
by a flow of fluid up the centre of the capillary. The free surface is relatively undeformed and appears to be 
close to its static shape, i.e. a spherical cap which meets the solid at the given contact angle. As one zooms 
in on the contact line region, it can be seen in Figure [6] that the domain is approximately a wedge of angle 
9 = 30°. It is also apparent that the FEM mesh has significant spatial resolution, so that, the dynamics of 
slip, on the (non-dimensional) length scale of /3^^ = 10^^, is well resolved as the smallest element has size 

kmn = O(10-«). 

For the same parameter values, in Figure [7] a comparison of the computed velocity Ug tangential to 
the liquid-solid and liquid-gas free surface with the local asymptotic result ( 43 1 is shown. Although the 
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Figure 7; Left: comparison of tlie velocity Us tangential to the free surface, curve 1, and solid surface, curve 2, with the local 
asymptotic result (dashed lines) from l |43| . The arclength along the interface from the contact line is s. Right: comparison 
of the normal stress along the liquid-solid interface A with the local asymptotic result (dashed line) from | |43[ l. The pressure 
constant po in l |43[ l has been chosen so that the numerical and asymptotic result coincide at the node adjacent to the contact 
line. 



computed free surface is curved, it is almost flat on the scale considered and hence Up, the radial velocity 
calculated in the asymptotic analysis, should agree well with Ug- In Figure [7| one can observe that the 
agreement between numerical and asymptotic results is good, with both velocities observed to be linear as 
the contact line is approached. Agreement between the computed and the asymptotic results is seen to 
occur on both surfaces for a distance from the contact line of s < ^S"^ = 10~^, i.e. inside the region of slip 
where the velocity starts to increase from its value at the contact line Ws = up to its far field value. 

The computed variable representing the normal stress A on the liquid-solid interface is seen in Figure [7] 
to agree very well with the asymptotic result (43); in particular, the normal stress is logarithmically singular 
as the contact line is approached. The agreement confirms that introducing this variable and incorporating 
it naturally into the weak form of the momentum equations has the desired effect. In particular, the level 
of approximation of this variable is sufficient. 

Thus, a comparison of our numerical results with local asymptotic analysis has shown that our code 
provides an accurate description of the fiow near the contact line. Next, we study the convergence of our 
scheme as the mesh is refined. 



5.3. Convergence of the new implementation 

We now consider whether the new implementation, in which the contact angle is imposed naturally into 
the weak form of the equations, ensures that 0c is equal to 9. In particular, high capillary number flow is 
considered as this is cited as the most problematic parameter regime in previous publications. 

In Figure|8j the influence that the capillary number has on the free surface shape is shown. For Ca < 10^^ 
the free surface shapes are graphically indistinguishable and represent the spherical cap one would expect 
a static meniscus to form. As one increases the capillary number, the apex of the free surface moves from 
below {za < 0) the contact line height to above it (z^ > 0). On the scale of the whole capillary radius 
the free surface does not seem heavily deformed; however, zooming in on the contact line region (Figure |8| 
reveals that the free surface actually makes an angle 9 = 30° in all of these cases. Therefore, at larger 
capillary numbers there is a small region of large curvature near the contact line. 

Each curve in Figure |8] was obtained using the new implementation and in each case the computed 
contact angle 9c was within 0.01° of 9 that was imposed. It appears that this can always be achieved when 
the mesh has sufficiently small elements near the contact line. To quantify this claim, we study the computed 
contact angle and, as a measure of the entire free surface shape, the difference in height H between the 
contact line z = and the apex z = Za for a range of meshes of increasing resolution, characterized by the 
smallest element size Imin, for Ca = 0.1. Note that the spatial resolution in the bulk of the mesh is the 
same in all of these simulations, it is only the resolution near the contact line which is altered. The results 
are presented in the table below. 
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Figure 8: Computed free surface shape on the scale of the capillary (left) and close to the contact line (right) as the capillary num- 
ber is varied with Re = 10, 6* = 30°, ^ = 10^ fixed. Curves 1-6 are for capillary numbers Ca = 10'^, lO"'', IQ-^, 10"^, 10-\ 1, 
respectively. 
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With increasing spatial resolution, the computed angle 9^ converges to the applied one 9 = 30° imposed 
via the weak formulation. For coarser meshes than those considered above, 9c exceeded 180° so that the 
free surface passes (unphysically) through the solid. It seems likely that such results are the reason that 
the contact angle is so often imposed as an essential condition, as in implementation (A-b). Critically, we 
see that the vertical distance H between the contact line and apex position is influenced by the resolution 
close to the contact line. So, errors generated by a poor approximation of the dynamics near the contact 
line greatly influence a measure of the entire free surface shape. 

We have shown that our results converge as the mesh size is reduced: using a specially designed mesh 
(i.e. high resolution with the appropriate increment of the element size) the computed contact angle 9c can 
be made arbitrarily close to the applied one 9, i.e. to the one featuring in the mathematical formulation of 
the problem. 

In the next section, we will examine the predictions of the numerical schemes previously used in the 
literature and will see that the conclusions from above also apply for implementation (A-a). What remains 
to be checked is whether by using implementation (A-b) , in which the contact angle is imposed as an essential 
condition, it is possible to obtain accurate results at a lower spatial resolution, i.e. at less computational 
cost, than considered above. This is the implicit assumption of numerical schemes which do not attempt to 
resolve the dynamics of slip and use (A-b) to ensure that the computed contact angle is always equal to the 
applied value. 

5.4- Examination of implementations (A-a) and (A-b) 

The same problem is now solved using the implementations previously considered in the literature. 
We found that implementation (A-a), with the one momentum equation applied at the contact line which 
is parallel to the solid, gives results which, for the accuracy considered in this paper, are indistinguishable 
from those obtained using the new implementation (B) . This means that for the particular flow configuration 
considered, the effect of reinstating the momentum equation normal to the solid surface is only to determine 
the function A, i.e. no additional accuracy is gained. Implementation (B) remains advantageous over (A-a), 
since there is no room for user 'input' at the contact line and the surfaces that are not aligned to coordinate 
axes require no additional work. It may be the case that the agreement between the results of the two 
implementations is due to the particular type of flow considered here with, for example, no component of 
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velocity normal to the solid. More complex flows in which there is a flux of mass into the solid, such as 
in the wetting of a porous medium, or in which the solid surface is curved may distinguish between the 
accuracy of the two methods. However, for the purposes of the present paper, the results from (A-a) and 
(B) may be presented together. 

Our focus is now on comparing the results from implementations (A-a) and (B) with those obtained 
using ( A-b) , where the contact angle is imposed as an essential condition, which is the most commonly used 
implementation in the literature. As shown in the table below, there is a significant difference in the results 
which these two methods produce. 
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By design, in implementation (A-b) the computed angle is equal to 9 for all meshes and, as hoped, 
this allows a solution to be found for coarser meshes. Furthermore, the value of H converges to the same 
value in all implementations, i.e. for a sufficiently refined mesh all implementations considered give the same 
result. However, imposing the contact angle in the weak form, as in approaches (A-a) and (B), is seen to 
create a formulation which converges to the correct solution significantly faster than approach (A-b). For 
example, for the third mesh, the percentage error in H using implementation (A-b) is more than a hundred 
times greater than that obtained from (A-a) and (B). This is despite the computed angle for (A-a) and (B) 
being more than 2° away from the imposed one, as shown in the first table, whereas (A-b) ensures equality 
between 6 and 9c- 

For the coarsest mesh considered, implementation (A-b) allows a solution to be obtained whose free 
surface is wildly different from the converged solution. This is despite the smallest element still being over 
one hundred times smaller than the radius of the capillary, i.e. being what one could easily interpret as 
'relatively small'. This serves as a representative example of the dangers of computing a solution on a mesh 
which is not specially designed to simulate dynamic wetting phenomena and then employing implementation 
(A-b) to 'compensate' for this and ensure a 'solution' is produced. 

As previously discussed, when the computed contact angle does not equal the one imposed in the weak 
form, it is tempting to impose the contact angle as an essential condition. This is because it is very easy 
to observe that the computed angle does not equal the applied one, more easy than, say, observing that 
the stress conditions on the free surface are not well satisfied near the contact line. However, what we have 
found is that the computed angle not equalling the imposed one should be interpreted by the user as a 
warning sign: the mesh is under-resolved. As seen from our results, attempting to bypass this problem by 
imposing the contact angle as an essential condition actually increases the global error of the scheme. 

Qualitatively, our conclusions are similar to the general remark made in the famous paper of Gresho & 
Lee [BH]i "Don't suppress the wiggles - They're telling you something", where it is warned that artificially 
suppressing spurious numerical artifacts merely shifts the generated error to another, less observable, part of 
the computed solution. We see, that when the contact angle is imposed naturally via the weak formulation, 
the deviation of 9c from 9 actually provides a simple quantitative error control on the numerical scheme. 
The error should not be suppressed by enforcing the angle as an essential condition, as this only exasperates 
errors in less easily observable parts of the scheme. One could insist that the computed angle is, say, within 
0.1° of the applied value and, from the results obtained above, in the parameter range considered, this 
would provide a very accurate numerical scheme. This is an idea which we now examine in order to provide 
practical quantitative recommendations for the spatial resolution required at given values of Ca and p. 

5.5. A practical guide to mesh design 

Using the aforementioned criterion that 9c is within a chosen distance from 9, we will now determine 
the required spatial resolution for given parameter values. At lower capillary numbers, roughly Ca < 10"^, 
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it is observed that the main constraint on a mesh-independent solution comes from the need to resolve the 
velocity field which varies on the (non-dimensional) length /?~^ near the contact line. The computed angle 
is seen to be within 0.1° of the applied value before this spatial resolution is achieved. 

At higher capillary numbers, where the interfacial curvature is high close to the contact line, the required 
spatial resolution must be determined. To quantify the results of the numerical experiments, we determine 
the size of the smallest element Imm that will ensure the computed contact angle is within 0.1° of the 
imposed one, in this case 9 — 30°. The results arc summarized in the table below. 







10" 


10^ 


10^ 


Ca 


10^2 


6 X 10-5 


6 X 10-*^ 


4 X 10-^ 


10-^ 


6 X 10-'^ 


5 X 10-' 


5 X 10-« 


1 
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7 X 10-« 
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A general trend can be seen that for higher Ca and /? it is required that l„iin is smaller, i.e. more 
spatial resolution is needed. Surprisingly, in this range of values, a simple empirical formula provides a 
good estimate for Imin- Imin = 5 X lO^^Ca^^ji^^ . It is important to remember that this is only applicable 
for large values of Ca, whereas for small Ca we still need to resolve the velocity field and this guarantees 
the closeness of Oc to 0. In practice, to ensure that the velocity field is resolved, we require that at least 
l-min < so that, finally, one needs 



I. 



5 X 10-^Ca-\l) 



(45) 



The linear dependence of Imin on Ca ^ can be rationalized by recalling that the local asymptotic 
result ( 44 ) predicted that the deviation of the free-surface from wedge-shaped (8 = 9) near the contact line 



is proportional to Ca/3. Thus, it is this product that determines the length scale on which the curvature 
of the free surface becomes significant and consequently on what scale the smallest element needs to be to 
ensure accurate results. We also considered the smallest element required for a 1° error in the computed 
contact angle and found that increasing the coefficient of proportionality by a factor of ten to 5 x 10-^ gave 
a good estimate for Imin- 

The above results can be used as a practical guide to the spatial resolution required to ensure accurate 
simulation of dynamic wetting flows. It is apparent that very small elements are required near the contact 
line which can only be obtained in a computationally tractable way by using a specially designed mesh. As 
we have seen, without this special attention to the contact line region the error in the computed solution 
can be huge. This will apply not only to the FEM but also to other numerical methods which, with a few 
exceptions (e.g. |30j ) . rarely incorporate the spatial resolution shown to be imperative for the recovery of 
meaningful results. 

In this section, we have fully resolved all issues surrounding the numerical implementation of a contact 
line at which an arbitrary angle is imposed. We now consider if the contact angle can be specified on the 
basis of physical considerations and how the resulting equations fit most naturally into the FEM framework 
which has been developed. 



6. FEM implementation of interfacial physics 

The contact angle is a boundary condition for the normal stress equation on the free surface, and, thus 
far, it has been considered as a prescribed input to the model. Now we will look into the 'extra' physics 
which actually determines its value, and examine the resulting effects on the flow. The contact angle is 
introduced into macroscopic fluid mechanics by the Young equation [69^ , which expresses the balance of the 
surface tension forces on the contact line acting tangentially to the solid surface: 

CTi cos 9d + (72^ 0. (46) 

Here, (Ti and (72 are the surface tensions acting on the contact line from the free surface and the liquid-solid 
interface, respectively, with the surface tension on the gas-solid interface being zero. Thus, the only way 
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in which the contact angle 9d can vary, i.e. be dynamic in the model and not as an artifact of numerics, 
is if at least one of the surface tensions varies from its equilibrium value die and cr2e, respectively. If the 
surface tensions become variables of the problem, this must be accounted for in the boundary conditions 
on the interfaces. Here, we allow for a dynamic surface tension and, being interested in the numerical 
implementation of the boundary conditions and for the sake of generality, do not concern ourselves with 
the mechanism for its variation (see [2] for such a mechanism). As it turns out, the equations on the free 
surface, and their implementation into the FEM scheme, already allow for gradients in surface tension. It is 
the equations on the liquid-solid interface where the influence of gradients in surface tension is not yet clear. 

We will now derive a generalized Navier condition by considering a balance of forces on an interface and 
an equation relating the torque acting on the interface to the velocity difference, i.e. the slip, across it. A 
balance of stresses acting on the liquid-solid interface is 

Ca n2 • (P+ - P-) + • - n^n^)) = 0, (47) 

where superscript -I- and — now denote, respectively, the value of a variable on the liquid-facing side of 
an interface, which is the value involved in the boundary to the Navier-Stokes equations (these variables 
have been referred to throughout the paper), and on the other side, in this case the solid- facing side of 
the interface]^ The stress on the solid-facing side of the liquid-solid interface is unknown and, if needed, is 
determined as part of the solution. Instead, the velocity on the solid-facing side of the liquid-solid interface 
U is known. A link between the two is obtained by relating the velocity difference across the interface to 
the difference in tangential stress acting on the interface, i.e. the torque on the interface, which gives 

\n2 ■ (P+ + P-) • (I - nana) - ^ (u+ - U,,) . (48) 

This condition is considered, for example, in [70] and later derived thermodynamically in [2]. Upon elimi- 



nation of the stress on the solid- facing side of the interface na • P from ( 48 1 using ( 47 ) , we arrive at the 
so-called generalized Navier condition p. 198], which has the form 

n2 • P+ • (I - nana) + .^Wa - /3 f u+ - Un) , (49) 



— V aa = ,.p| -^iiy, 

and can see that, when the surface tension is a constant, this generalized condition reduces to the standard 
Navier condition 



Using (491, we now have the contribution to the momentum residuals (21) on the liquid-solid interface 
given by 

A (e^ • n) /3 • (u|| - U||) - i^r^^a ■ Wa 
This new formulation has been implemented into our numerical scheme 



A2 



dA2. (50) 



6.1. Verification of the method 

The numerical scheme containing the additional physics has been thoroughly tested and gives excellent 
results. To demonstrate this, let us examine the same axisymmetric flow in a capillary as considered in |5.2[ 
taking Ca = 0.1, Re — 10, P = 10^ and prescribing the surface tension on the interfaces to be 

0-1 = (Tie = 1, era = (T2e + ^ exp(/32;). (51) 

Here, CTae — — -\/3/2 is taken to be the value of the equilibrium surface tension on the liquid-solid interface, 
which is the asymptotic limit of ua in the far field (as z — ^ — oo) and ensures that in equilibrium the Young 



^The same equation applie s on the liquid-gas interface which, when combined with a condition on the gas-facing side of the 
interface n ■ =0, yielded ISSl. 
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equation (46) gives that cos = 
line z = 0, from ( 51 ) one has (T2 



—<J2e/<^ie = \/3/2, SO that 9e = 30°. In the dynamic case, at the contact 
(1 — \/3)/2 so that the dynamic contact angle 9^ = arccos(— (T2) = 68.53°. 
Although the formulated problem contains artificially devised surface tension distributions, it is still of 
interest to see the influence which these have on the flow field. In Figure |9j we show streamlines for this 
flow in the region where both the tangential stress and the surface tension gradients contribute to slip at 
the liquid-solid interface. It is seen that the surface tension gradients can be strong enough to reverse the 
flow close to the contact line. 



For axisymmetric flow, equation ( 49 ) can be written down in scalar form as 



dw 1 d(T2 

dr 2Ca dz 



/3(w; + l), 



(52) 



and in Figure [9] we plot the behaviour of the individual terms from this equation. The equation is well 
satisfled with curve 1 and curve 2 adding to give curve 3. From this figure the relative strengths of the 
factors which cause slip can be seen: it is the surface tension gradient which is largest for s < 10^^ and 
hence it is able to drive a flow reversal. 

In summary, it has been shown that our new framework is easily generalizable and, in particular, allows 
interfacial physics associated with dynamic wetting to be naturally incorporated into it. 



7. Discussion 

As has been discussed, a number of different FEM implementations have been used to simulate dynamic 
wetting phenomena. We have shown that to guarantee accurate results, particularly at higher capillary 
numbers, sufficient spatial resolution of the computational mesh must be used near the contact line. When 
this is the case, many of the seemingly different implementations, including our own, converge to give the 
same result. Despite these similarities, approaches in which the weak formulation is used to enforce the 
contact angle are seen to converge to the correct solution faster. An apparent 'window' for a numerical error 
imposing the contact angle using the weak formulation is seen to actually be an advantage: the difference 
between the applied angle and the computed one provides a simple quantitative test on the accuracy of the 
scheme. This constraint is lost if the angle is imposed as an essential condition whilst the error is shifted 
elsewhere in the code but not removed. In particular, imposing the angle as an essential condition when 
using a mesh which does not incorporate sufficient extra spatial resolution near the contact line can lead to 
huge global errors in the free surface shape. 
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In the new implementation, the momentum equations are applied everywhere and any surface which is 
curved or does not align with coordinate axes causes no additional complexity. A key advantage of our 
framework is that, since there are no 'loose ends' at the contact line, there are no choices left for a user in 
how to implement the contact line region. This, together with the obtained quantitative criterion for the 
mesh, ensures that the implementation of dynamic wetting phenomena into the FEM satisfies the mantra 
in [711 p. 35] that, once the weak form is derived and the choice of mesh design, element type, etc are taken, 
"the rest of the recipe is well defined - just turn the crank to generate the complete spatial approximation". 

It has been shown that when the weak form is used to impose the contact angle, the computed value 
on an under-resolved mesh can vary significantly as the mesh resolution is altered near the contact line, 
and to interpret the computed value, a mesh-dependent quantity, as the dynamic contact angle would be 
fundamentally incorrect. If one wishes to define an 'apparent' angle, that is some measure of the free-surface 
shape away from the contact line, then this should be a well-defined variable, see e.g. which converges 
to a specific value as the mesh is refined. In phenomena where the dynamic contact angle starts from an 
out-of-equilibrium value, e.g. in drop impact and spreading phenomena, where the angle starts at 180°, it is 
a fundamental modelling mistake to prescribe the contact angle to be equal to its equilibrium value. In this 
case, the applied angle cannot be reconciled with the initial free-surface shape in principle, and the code 
will be unable to provide mesh-independent solutions for early time. 

One may hope that the restrictions on mesh design described in this paper could be relaxed by removing 
the region near the contact line and, if an appropriate limit is available, incorporating some local asymptotic 
results there. This approach has been considered, for example, in [7^ [75| and is discussed in [73]. The idea 
is to incorporate the small capillary number asymptotic results of [75] . which relate the contact angle to 
the free surface shape at a given distance from the contact line. In both investigations, it is seen that such 
an approach fails for Ca > 0.1. Therefore, there is a limited range of validity to these numerical schemes 
which prevents them from being used to simulate flows where the capillary number is of order one or higher 
regularly encountered in processes of technological interest, e.g. the impact and spreading of liquid drops 
ejected from ink-jet printers. Furthermore, the approach taken in these articles relies on the assumption 
that the flow field has no influence on the dynamic contact angle. This has been shown by experiments 
to be incorrect, e.g. in [TH [TH]. The model aimed at describing this effect has to be based on the Young 
equation that actually introduces the contact angle as a macroscopic fluid-mechanical concept and hence its 
numerical implementation has to rely on incorporating gradients of the surface tension along the contacting 
interfaces [2]. We have shown that the developed FEM framework allows such a generalization to be made 
in a naturally way without any deviation from the FEM methodology. 



8. Appendix. FEM implementation for 2D and 3D axisymmetric Q.ow 

The Appendix provides a user-friendly ready-to-use guide to the specific implementation of our FEM 
framework for a dynamic wetting process. This guide enables one to easily reproduce all the results presented 
in the paper. For a more detailed exposition about using the FEM to solve fiuid flows the reader is referred 
to j34j[71j, whilst specific information on free surface fiows, including alternative mesh designs, can be found 
in [2Tti23| [35 l [76] . For a beginner's guide to developing and structuring finite element code we suggest [77] . 

Two-dimensional and three-dimensional axisymmetric flows are considered simultaneously by using a 
variable n which takes the values n = in the former case and n = 1 in the latter. Two-dimensional flow 
occurs in a Cartesian coordinate system (r, z) whilst axisymmetric fiow is in a cylindrical polar coordinate 
system (r, z, i?), where r is now the radial coordinate and i? is the azimuthal angular coordinate around the 
z-axis about which symmetry is assumed (Figure |4]). For both coordinate systems, the governing equations 
are solved in the (r, z)-plane. 

The Appendix is organized as follows. First, in Section [8. 1[ the tensorial expressions in the main body of 
the text are specified for the coordinate systems considered. Section [8?2] describes the arbitrary Lagrangian- 
Eulerian mesh design, based on the method of spines, which tessellates the domain into elements whose 



position becomes dependent upon the free surface shape. In Section 8.3 , the interpolating functions are 
constructed piecewise over each clement so that both the functions and the residuals are calculated element- 
by-element. Expressions for the residuals on an arbitrary element are determined. Section |8.4| shows that 
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by mapping each element in the computational domain onto a master element the integral expressions 
from Section 8.3 can be calculated in a systematic way. In Section 8.5 a method for solving the resulting 
set of algebraic equations using a Newton-Raphson method is outlined. Specific attention is given to the 
construction of the Jacobian which is required in the iteration process. Finally, in Section |8.6[ additional 
details of the numerical scheme are given. 

The Appendix provides information of our own specific implementation of the framework described in 
the main body of the paper. Many alternatives exist (e.g. different element types) which may produce 
equally accurate results but here we present what we have used and tested, which for a 'practitioner' gives 
a ready-to-use algorithm. 

8.1. Scalar expressions 

We will first give coordinate-specific expressions for both two-dimensional flow (n = 0) and three- 
dimensional axisymmetric flow {n = 1) in the integral form and then consider the differential terms in these 
integrals. 

8.1.1. Integral expressions 

The integral of a function / over a volume V , corresponding to an appropriate portion of the (r, z)-plane, 

is 

'' f dV ^ ( I f (27rr)" dr dz, (53) 



V Jz Jr 

where integration over the d coordinate for the axisymmetric case has yielded the (27rr)" term. Surfaces are 
parameterized in terms of the arclength s = s in the (r, z)-plane and, in the axisymmetric case, the angular 
coordinate t = d. The integral of a function / over the surface A, which is a line in the (r, z)-plane, will 
then be 

^ f dA^ f f (2^r)" ds, (54) 



for appropriate values of s. Finally, integration over the contact line contour Cd, which in the (r, z)-plane 
is a point (rc, Zc), will give 

f dA={2nr,rf{r,,z,). (55) 



Note, the factor (27r)"' will appear from all integrals and is thus cancelled throughout the equations given 
in Section [831 



8.1.2. Bulk expressions 

Using unit basis vectors of the coordinate system and e^, we decompose the bulk velocity into 
components along the coordinate axes as u = uer + wez and the components of the solid substrate speed 
as U = Uej. + WGz- We can now derive the coordinate-specific expressions for the differential terms that 
form the integrands. 

The bulk gradient operator is 

9 n d d 
V e^— +ei,-— +6^ — . (56) 
or r ov oz 

Then, the gradient of a vector A = ar(r, z)er + az{r, z)e;,, which will be a second order tensor, will be 
required and is given by 

da-r , da-r , nOr ^ da^ da^ 

VA = ^— e^e^ + ^— e^e^ H e^^e^, -f — — e^e^ + — — e^e^, (57) 

or oz r or oz 

whilst the divergence of the vector A, which is the trace of its gradient, is 

Or r Oz 
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The stress tensor is given by 



P = Pr^erCr + Prz^rez + P^^B^B^ + Ezre^Br + PzzBzBz, (59) 

with components 

„ r,9u „ „ dw du „ 2nu „ ^dw 

Prr = -p+ 2 — , Prz = Pzr ^ ^ + ^ = -p + , Pzz = -p + 2 — . 60 

or or oz r oz 



Using (57) with A — (piBr and then A = (t)iBz combined with (59), an expression for V(0ieQ) : P, which 
is required in the bulk momentum residuals, is obtained: 

V{<i)iBr) :P = ^Prr + ^Prz ~i P-di) , V (0,6^ ) : P = — P^^ + ^ -P^^ • (61) 

or oz r or oz 

8.1.3. Boundary expressions 

In the (r, z)-plane, surfaces are parameterized in terms of the arclength s, which increases away from the 
contact line, so that the coordinates of a point on a surface are {r{s), z{s)) or in the case of the axisymmetric 
system (r(s), z(s), i?). The tangent vector in the (r, z)-plane, which points in the direction of increasing s, 
is t = trB^ + tzBz and the vector normal to the surface, which points into the liquid, is n = n.rBj. + UzBz. 

The vectors normal and tangent to the surface are then 

(^^^^^)= (^j!|!;l)i/2 ^ {nr,nz) = ±{-tz,tr), (62) 

where the prime denotes partial differentiation with respect to s and the sign is chosen so that the normal 
vector will point into the liquid. In axisymmetric coordinates, the second tangent vector, perpendicular to 
the one in the (r, z)-plane, is the unit basis vector in the azimuthal direction e^. 

Now, the tensor (I — nn), which is used to extract the tangential component of vectors, has the form 

(I — nn) = tt + rte^e^, (63) 

and hence we can also find the surface gradient operator 

= (I - nn) • V = t (t • V) + ne,, (e^ . V) = t— + ^ — . (64) 

OS r ox) 

Therefore, in particular, expressions for the terms associated with the stress on the free surface are 

V • ((/-^ej = (t • e^) — H , (65) 

OS r 

V^.(0,ej = (t.e,)^, (66) 

OS 

QBj. 9Bz 

where we have used -—— — b§, = 0. 

ov Ov 

All of the components which appear in the residual equations have now been derived. Next we consider 
how the finite element procedure begins by tessellating the domain into elements. 



8.2. Method of spines mesh design 

We use a well known arbitrary Lagrangian-Eulerian approach, known as the method of spines, which 
has been successfully applied to many free surface flows, e.g. [551 US]- This method, which is described in 
detail in [3T], builds on the boundary location method proposed in [33]. In the method, nodes which define 
the free surface are located at the end of a line on which bulk nodes are attached. These are the so called 
spines which form the skeleton of the mesh and usually run between a solid surface and the free surface, 
with nodes often spaced equally along each spine. As an iteration proceeds, all the free surface nodes evolve 
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to new positions and hence the spine, which is attached to the free surface node, moves and drags the bulk 
nodes accordingly. Elements are tesselated around the nodes and therefore their position is dependent on 
the free surface shape. 

In the (r, z) computational domain the free surface is a line parameterized by the arclength s. Then, 
the position of the free surface can be defined by a function of one variable h = h{s) which depends on the 
mesh design chosen. Finding h at each node hi on the free surface (i = 1-iVi) gives the free surface shape. 
An example of the free-surface parameterization is given in §3.4[ and more specific schemes are presented 
below. 

The ideal spines to capture the dynamics near the contact line, where the geometry is wedge-like, are 
arcs of circles whose centre is the contact line (Figure [5]) and whose arc angles are the /I's. We will show 
that this can be easily achieved by introducing a polar coordinate system centred at the contact line. The 
drawback of this polar method is that it can be difficult to match the mesh design near the contact line 
with the geometry of the computational domain, which can often involve straight boundaries, further away. 
To overcome this problem, we also present a method based on the bipolar coordinate system which creates 
circular spines near the contact line and straight spines further away 

To ensure the computational problem is tractable and that the smallest scales in the model are well 
resolved, a graded mesh is used with the smallest element near the contact line. This is achieved by defining 
each spine to be a given distance from the contact line. In particular, if we choose to have Nk spines, with 
fc = 1 corresponding to the contact-line point and k — 2 the first proper spine, i.e. an arc, and increase the 
distance between adjacent spines by a ratio r, then the distance Rk of spine k from the contact line is 

k— 1 2^ 

Rk - Rmax I ^ = ^-^k , (67) 

where Rmax is the distance of the final spine (fc = Nk) from the contact line. In this paper, we used r = 1.07, 



with Nk varied to create meshes of differing spatial resolutions. Expression (67) is used for each of the two 
methods presented. For simplicity, when discussing the mesh design we shall consider, as is the case for the 
capillary, a free surface meeting a fiat solid (r = 1) at a contact line (1, Zc). We note that in the actual code 
the contact line is at z = Zc and is free to move. Only once we have obtained a converged solution do we 
map z — > z — Zc so that the contact line is at (1,0) as in the presented results. 

8.2.1. Polar spines 

A polar coordinate system (i?, ip) is created around the contact line, which is related to the (r, z) coor- 
dinates by 

(r, z) = (1 — i?siniy9, z^ — Rcosip) . 



Therefore, given the polar coordinates of a node, we can obtain its actual position in the computational 
domain. The free surface unknown at the contact line = 1 is the contact line's height hi — z^. Spines 
k = 2-Nk are arcs of a circle of radius Rk and arc angle (fk from the solid surface to the free surface. The 
angle of the arc is the free surface unknown which is to be determined as part of the solution, i.e. hk — fk- 
Nodes are then placed uniformly across each spine, so that node m out of M nodes on a given spine k is at 
angle 

m — 1 



fni,k = Vk [jl^ j ' (69) 

and the subscript m, k refers to the fact that this is the m-th node of spine k. Then, for given to, fc a node 
in the domain is dependent upon the vertical coordinate of the contact line hi = Zc, which determines the 



(known) distance of spine k from the contact line calculated from (671, and the angular coordinate hk = fk 



which determines how nodes are distributed across that arc (Figure 10 1 



8.2.2. Bipolar spines 

Far away from the contact line, in the region of the capillary between the solid surface and the axis of 
symmetry, we would like to use straight spines to match with the capillary's geometry. Such spines cannot be 
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matched to circular spines emanating from the contact hne region. However, by using the bipolar coordinate 
system, the spines near the contact line can remain circular whilst those further away, at a distance chosen 
by the user, are straight. To achieve this, the bipolar coordinates are first related to arbitrary Cartesian 
coordinates. Then, an appropriate mapping is used to embed the mesh in the computational domain. 
The bipolar coordinate system (x, C) is related to Cartesian coordinates (x, y) by 

r sinhx , sinC , . 

^ = / — u — , — 7' y = / — u — , — 7' 

cosh X + cos C cosh x + cos C 

where {x, y) = (/, 0) is the focus of the bipolar coordinate system. The bipolar coordinate x is used to define 
spines, as R was previously used, with spine k having x — Xk and being a circle with radius and centre 
given, respectively, by 



smhxfc Vtanhxfc / 

so that as Xfe — ^ oo the spines' centres approach the focus and their radius tends to zero whilst as Xfc — > 
the radius of the circles tends to infinity, i.e. the radius of curvature tends to zero, and we approach straight 



lines at j: = (Figure 11 ). 



We want to map the bipolar coordinate system in {x, y) with focus at (/, 0) and straight spines at a; = 
onto the computational domain (r, z) with the focus at the contact line (1, Zc), where spines of a small radius 
are required, and straight spines at a chosen distance Rmax from the contact line so that / = Rmax- Then, 
on a given spine k with x = Xfcj varying C, will take us from the (fixed) solid surface at ^ = ^s, with to 
be determined, up to the free surface at C = Cfci i-^- C takes the role which Lp previously played for the polar 
method. Notably, here we have not assumed that the solid surface occurs at y = which would correspond 
to C = 0. Choosing Q. to be the angle at which the a;-axis meets the r-coordinate this can be achieved by 
the transformation 

( r~ro\_f cos{n) - sm{n) \ f x \ , . 

\z-z^ )-\ sin(f^) cos(f^) )\y )' ^ ' 

where the origin of the {x,y) coordinate system is placed at (rpj^o) = (1 ~ Rmax cosft, Zc — RmaxSinfl). 



Therefore, given values of (XiC): '^c can determine {x,y) from (70) and hence a position in the (r, z)-plane 



from ( 72 ) . One can still use mjj to define the distance along the a:-axis {( = 0) of a given spine from the 



contact line, except one must now use (70) and (72) to determine the appropriate value of Xk from Rk- For 



our particular geometry we take Rmax such that the y-axis passes through the apex of the free surface at 
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Figure 11: Left; bipolar spines in a Cartesian frame highlighting the coordinates at the start and end of the fifth spine. Right: 
embedding of the bipolar mesh system into the computational domain. 



z = Za and choose the angle i7 such that the distance along y to the apex is the same as the distance along 
y to the solid (Figure 111. 

Outside the bipolar region, one is left with a straight-sided domain whose upper boundary is dependent 
on the position of the final, straight, bipolar spine. The spines in this region are chosen to run from the solid 
surface to the axis of symmetry, both of which are fixed, so their positions may be easily defined in terms 
of the final bipolar spine. What has been described is only one possible mesh, which is shown in Figure [Sj 
that can be designed by mapping the bipolar coordinate system onto the computational domain and is seen 
to work well for the range of parameters considered in this paper. 



8.2.3. Remarks on mesh design 

The aforementioned procedures create elements throughout the domain which have curved sides. How- 
ever, only the element sides which form the free surface need to be curved and so, if required, to straighten 
all the other element sides one can loop through every node and define every midpoint node to lie exactly 
half way between its two adjacent vertex nodes. The drawback of this approach is that the midpoint nodes 
become dependent on the position of the adjacent spines. This means that the nodal positions become 
dependent on, at most: the position of the contact line; the apex position, which is used to define in (72); 
and on the free surface unknowns of the adjacent spines. 

The polar and the bipolar meshes can be used local to the contact line with a different mesh design, or 
even a different numerical method [78], taking over further away. In our case, we used straight spines for the 
rest of the capillary as the geometry is simple; however, it would be possible to use unstructured meshing 
away from the contact line which would allow a greater degree of flexibility. 



8.3. Finite element procedure 

Having derived the expressions which will be required in the residuals and described a particular mesh 
design, we proceed to outline the finite element method. The method is not restricted to the capillary 
geometry and our only assumption about the free surface is that in the (r, 2:)-plane, i.e. the computational 
domain, it is a line parameterized by arclength s with position defined by a function of one variable h ~ h{s) 
which depends on the particular mesh design. In other words, to determine the free surface shape we must 
determine the value of hi at every node on the free surface i = l-7Vi . 

The domain is tesselated into e = l-TVg non-overlapping subdomains of area Ee, the elements, which 
in our case will be two-dimensional curved-sided triangles whose positions are defined by the spines of the 
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Figure 12: Two elements adjacent to the contact line showing the relationship between the element numbers e = 1, 2, the global 
nodes (which are given on the outside of the elements) and the local node numbers (on the inside of the elements) . Also shown 
are the first line element on the free surface si^ which is a side of element e = 2 and the first line element on the solid surface 
S2e which is a side of element e = 1. 



mesh and hence by the values of free surface unknowns hi. The boundary of the computational domain is 
composed of one-dimensional elements of arclength Se formed from the sides of the bulk elements which are 
adjacent to the boundary. We refer to element level quantities as local and those defined over the entire 
domain, which were used in ^jS] as global. Each element contains a set number of local nodes and it is the 
set of all local nodes which form the global nodes referred to in f|3] We have used Roman letters (i,j) for 
global quantities and italicized letters {i,j,k,l) for local ones. In the FEM, one must store a function I 
which relates each local node i in each element e to its global node number i so that i = /(e, i). Each global 



node often belongs to more than one element, see Figure 12 where, for example 7 — /(1, 3) — I{2, 1). 

In this section, it will be shown how the global functions and residuals are constructed in a piecewise 
manner from local quantities, by ensuring that, as required in the FEM, the interpolating function associated 
with a given node is constructed to be zero outside the elements to which that node belongs. This allows us 
to derive expressions for the residuals in an arbitrary element which, after summing up contributions from 
every element in the domain, will form a set of algebraic equations whose coefficients are integrals to be 
calculated using the methods described in Section |8.4[ 

8.3.1. Construction of the interpolating functions 

The global interpolating functions are constructed piecewise from local interpolating functions which are 
defined across elements as low order polynomial functions. A global interpolating function, say is zero 
in all elements which do not contain global node i. In an element e whose local node i is global node i, the 
global function is equal to the local one (pi. For example, the global velocity interpolating function in an 
arbitrary element e is: 

otherwise. ^''^I 
To avoid having to construct each local interpolating function in an element-specific manner, they are 



constructed on a master element with coordinates {£,,ri), see Figure 13 so that (f>i = (t>i{S,,r])^ tpi = ipi^S.yrj). 
Mixed interpolation is achieved, to ensure the Ladyzhenskaya-Babuska-Brezzi [40| condition is satisfied, by 
using the V6P3 Taylor-Hood triangular element which approximates velocity by means of bi-quadratic local 
interpolating functions ipiHi v) (* — and pressure using bi-linear ones v) = 1-3). As discussed, 
the interpolating functions are constructed to take the value one at local node i and zero at all others 
associated with that function, so that 

^i = ^^-^> ^2 = -^^, V'a = (74) 
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Figure 13: The V6P3 master element with local coordinates (5, rj) which is mapped into the physical (r, z) space via the 
isoparametric mapping l |76| . The numbered circles represent nodes at which the velocity and pressure are to be found whilst 
the squares are velocity only nodes. 



(75) 



and 

= ^i(27/;i-l), 02 = 1/^2(2^2-1), 03 = M^^3-l), 

04 = 47/;iV'3, 05 = 40>2V'l, 06 = 4V'3'02- 

In a given element e, whose nodes in the global domain are located at (rj, Zj), (j = 1-6), the relationship 
between a position in the global domain (r, z) and a position in the master element rj) is given by 



6 



6 

E 



Zj(/l)0j(^,77), 



(76) 



where the dependence of all of the nodal positions on h, i.e. the free surface unknowns, is shown here 
explicitly but henceforth assumed. Being able to find the global coordinates from only the nodal positions 
and the interpolating functions gives the mapping the property of being isoparametric. 

Next, we form the one-dimensional local interpolating functions which are required for the approximation 
of boundary terms. For example, on the side of the master element (Figure 131 with nodes j = 2,6, 3 a 
one-dimensional master element rj = —1 and ^ G [—1,1] can be defined on which a one-dimensional quadratic 
(as opposed to bi-quadratic) interpolating function 0ij(^) = 4'jiS,iV ~ —1) is obtained from (75): 

1 + e 



MO 



^^1,3(0 



'^1,2(0 = V'1,2(2V'1.2 - 1), 01,3(0 = V'1,3(2V'1,3 - 1), 01,6(0 = 4V'1,3V'1,2. 

Then, the coordinates in the global domain on this side of the element are 



ri = 



3=2.6,3 



r3<PiA0 ^1=2. ^3M0- (77) 

1 J=2,6,3 

By ensuring that elements whose sides form the free surface always have local nodes 2,6,3 associated with 



them, we can use (771 to define the free surface shape. Similarly, by ensuring that the local nodes on the 
solid surface are always j = 1, 5, 2, we can define (t>2.j(j]) = 0i(C = — 1, v)- 

Having constructed the required interpolating functions, over an element e the following functional forms 
are obtained 



{u,w) = ^{uj,Wj)(j)j{^,'q), 
j=i 



J=i 
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(r,z) 



6 

i=i 



r3^Zj)(l)j{^,ri). 



(78) 



The surface tension is also approximated on the surfaces quadratically as 



3=2,6,3 j = l,5,2 

Note that, in this article the nodal values of the surface tension were prescribed. 

The global residuals (the i?i's) presented in ^ which involve integrals over the entire domain, are formed 
by summing up local residuals (the i?e,i's), denoted with a superscript e to indicate which element the local 
residual is calculated in, obtained by integrating over each of the elements in the global domain. Let 
Ni be the set of numbers of those elements that form part of the free surface and N2 be the numbers of 
elements forming a part of the liquid-solid interface. Then, the continuity of mass, momentum equations, 
kinematic equation on the free surface and impermeability equation at the liquid-solid interface are: 

Ne 6 We 3 Afe A^e 

e=l i=l e=l i=l e=l i=2,6,3 e=l i=l,5,2 

i=-f(e:0 i=I{e,i) i=/(e,i), eeJVi i=I{e,i), eSiVa 

where the constraint under the summation symbols ensures that, through i — I{e,i), the local residuals are 
assigned to the correct global ones and, via e G N, that the surface equations are only evaluated when an 
element is adjacent to the boundary of the domain. 



3.3 



The local residuals in (80) are simply the equations of ^3.3 consider over a subdomain of V. Notably, 



it was shown that the manipulation of the momentum residuals resulted in terms which 

represented the stress acting on the bounding surface of the domain. One may expect that these terms 
will contribute to the momentum residuals on the boundaries of each element. However, at nodes on inter- 



element boundaries, such as global node 3 in Figure 12 there will be a boundary contribution to the global 



residual from two adjacent elements, with normal vectors of opposite signs, both representing the stress on 
the edge of the element. By cancelling these terms at inter-element boundaries continuity of stress is ensured 
across the interior of the domain. Then, surface contributions to the local residuals only occur when one of 
the sides of an element e forms part of the free surface e G Ni ot the solid surface e G N2. In the same way, 
by doing nothing, we can ensure continuity of surface stress across all free surface elements: this actually 
amounts to continuity of the tangent when the surface tension is continuous. Contributions from the edge 
of a free surface element only occur at the contact line. 

Taking the equations of ^|3]over an arbitrary element e, with area and boundary Sg, and using the 



expressions in Section 8.1 we now derive the local residuals required in (80 1. The expressions appearing in 



these equations, which are usually integrals of products of the interpolating functions and their derivatives 



over the element or its boundary, are given in Section 8.3.3 



8.3.2. Element-level residuals 

For an arbitrary element e, with local nodes i = 1-6, contributions from the bulk of the element E^, to 
the momentum residuals are 

R^f = {Re A,, [u, w,h)+ Klj (h)) + Klf {h)w, + C}^ [h)pu j = 1-6 fc = 1-3, (81) 
=Kfl {h)uj + {Re A,, [u, w,h) + K^J (h)) Wj + Cf^ {h)pk j = 1-6 k = 1-3, (82) 

where summation over repeated indices is henceforth assumed. The K terms are associated with viscous 
forces, the A terms are from the nonlinear convective terms and the C terms are with pressure forces. 
Expressions inside brackets indicate the dependencies of each term on other functions in the problem, i.e. 
they show the nonlinearity of the problem. 

For i = 1-3 we have the incompressibility residuals 

R% - C], {h)uj + Cl (K)w, J = 1-6. (83) 
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If an element e forms a part of the free surface, e € Ni, then for i ~ 2,6, 3, i.e. the free surface nodes, 
there are additional contributions to the momentum equations from capillary stress terms so that 

<r=<r + Fl,ih)a,, J = 2, 6, 3, (84) 



=Kf + 4 (M^i.. 3 = 2, 6, 3. (85) 



I] \ 

Additionally, for i = 2, 6, 3 there is kinematic equation of impermeability 

i?5 = 4 {h)u, + 4 {h)w, J = 2, 6, 3. (86) 

If an element e forms a part of the liquid-solid interface, e G N2, then for i = 1, 5, 2 there are additional 
contributions to the momentum equations of 



R^T + NlHh)^, + NlKh)w, + L\h)\, + SUh)a2,, 3 = 1, 5, 2, (87) 



<r + N}l{h)u, + ^'(/i)«, + i'(/i)A, + 4(/i)f72,, J = 1, 5, 2, 

where the terms come from either the Navier condition Q or its generalization ( 49 1 , whilst the S terms are 
associated with gradients in surface tension which are only accounted for in the generalized Navier condition 
(49 1. The L terms are associated with the normal stress on the interface. Additionally, for i = 1,5,2 one 
has the impermeability equation ([6]) in the form 

= 4 ('^)%- + 4 (^)^^ J = 1, 5, 2. (89) 

If an element e forms the part of the free surface e G iVi which is adjacent to the contact line, so that 
local node i = 2 is the contact line node, then there is an additional contribution 

Rt^ (90) 
R't^i' =Rti' +T\h,Q), (91) 

which is where the contact angle is imposed into the weak formulation. 

8.3.3. Required matrices 

The viscous terms, see Section [8?Tj are 



where 

l^'-^r-dE^, "i^r-dE^, (93) 

Je^ dyk dyi ^ Je^ r2 

with Hi = r, 1/2 = z. Nonlinear convective terms are 

Aij{u, w) = aijkiUk + aijk2Wk k = 1-6, (94) 

where 

a^Jkl = [ <|>^c|)k^ r" dE,. (95) 
The C matrices represent pressure terms whilst their transpose is required in the continuity of mass equation: 



Kll = 2fc,jii + hj22 + 2n fc[^., Kl^ = %2i, Kfl = fc„-i2, Kf^ = h^n + 2fc„-22, (92) 
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On a free surface (one-dimensional) element sie, the capillary stress terms are 



Ca 



ds 



r 



1 



ds 



(97) 



The components of the vectors tangential t — (1^,1^) and normal n — {rir, Uz) to a surface have been derived 
in Section [01 

On a (one-dimensional) element S2e on the liquid-solid interface, the tangential stress terms from the 
Navier-slip boundary condition are 



N, 



kl 



P / 0: 



t2kt2l ds2e, 



(98) 



with t2i = t2r and i22 = t2z- In the generalized Navier condition (49 1, additional terms associated with 
gradients in surface tension are 



Sh 



1 



2Ca 



'2j 



ds 



t2r r" dS2e, 



and the normal stress terms give 

Llj = / 4>2,i(t>2,jn2r dS2e 



St 



- 



1 



2Ca 



»2j 



ds 



t2z r" dS2e 



t)2,i(P2,]n2z r 



'' ds 



2e- 



(99) 



(100) 



At the contact line (rc, z^) contributions will occur at local node i = 2 of the first free surface element 
so that 



[{t2r cose + n2r sin 0)r^], 



[{t2z cos9 + n2z sin 6*) r"] , 



(101) 



Ca ^' Ca 
where 6 is either prescribed or determined from the Young equation 9 = arccos(— CT2,i=2/o'i.i=2), with the 
surface tension values evaluated at the contact line. 

On both the free surface (7 — 1) and the liquid-solid interface (7 = 2), expressions for the impermeability 
equation in surface element s-^g are 



If. 



(102) 



We have now derived all the expressions required. In the next section we present a systematic way 
of determining these integrals and hence forming a set of nonlinear algebraic equations with numerical 
coefficients. 



8.4- Mapping local integrals onto the master element for evaluation 

We have a set of algebraic equations, that is the momentum equations (81 )-(82 1 with contributions from 
the free surface (84)-(85), solid surface (87)-(88l and contact line (90|-(91); the continuity equation (83); the 
kinematic condition on the free surface (86 1; and the impermeability condition on the solid surface (89). The 
coefficients of the functions' nodal unknowns in these equations are the integrals over deformed elements 
in the global domain given in Section 8.3.3 These integrals are most easily determined in a systematic 
way by mapping them onto the master element (Figure 13 1, over which we already have constructed the 
interpolating functions, for evaluation. Therefore, each integral will be over the same geometrically simple 
domain and, as described in Section |8.4.4[ quadrature methods may be used to accurately approximate 
them. 

From the mapping (76) we have a relationship of the form r — r(^, rj), z — z(^, 77) which will be required 
to transform our integrals. Derivatives of the global coordinates with respect to the master element's 
coordinates will also be required and are easily expressed as 

- .a n03^ 

di ~ ' dC dTi~ ' 977 ' di~ ' di' dv~ ' dr^ ^~ ^ ' 

Therefore, we would like all integrands to contain the interpolating functions and their derivatives with 
respect to the master element's coordinates, as these are quantities are easily calculable. 



34 



8.4-1- Evaluating local 'volume' integrals 

The integral of a function / over a (two-dimensional) element in the computational domain is trans- 
formed into an integral over the master element as follows 



/ dE„ = 

where the Jacobian of the isoparametric mapping is 



=1 


=-11 




-1 


I is 




/ dr 


dr 
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drj 

dz 




dr] 



f det Je dr], 



(104) 



so that, using (1031, its determinant is 



det Je — 



dr dz dr dz 



d^ dr] dr] d£, 
where we have introduced the useful matrix 



Eij r^ Zj 



(105) 



i,.] = 1-6, 



(106) 



dS, dr] 



dr] dS, 



(107) 



We can ensure that the determinant is positive, and therefore do not require its modulus, by making sure 
that the local numbering of our master element and of our elements in the global domain are both, say. 



anti-clockwise (Figure 13). 

The derivatives of the master element's coordinates with respect to the global ones may be found by 
noting that the inverse of the Jacobian may be obtained either by considering ^ = ^(r, z), r] = r](r, z) or by 



directly inverting (1051, so that 
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d^ 



(108) 



d^ d z 
Then, for example, — = — / det 



dr dr] ' 

To determine the derivatives of a variable with respect to global coordinates, in terms of the (known) 



element's global coordinates and (known) derivatives of the local interpolating functions (74) and (751 with 



respect to local coordinates, we use the chain rule and the expressions in ( 108 1 to obtain 

TikZk 



dr 



dz 



^ d^d^ 
d^ dr 
d^d^ 
dC dz 



d<j}i dr] 
dr) dr 
d(l)i dr] _ 
dr] dz 



det Je 
Tjkrk 
det Je 



k = 1-6, 
k = 1-6. 



(109) 
(110) 



We have now obtained all the required expressions to allow us to evaluate integrals over a standard 
master element instead of over the deformed global elements. Next, we consider the same procedure for the 
'surface' integrals. 



8-4-2- Evaluating local 'surface' integrals 

In the computational domain, the 'surface' is divided into line elements Se which are the curved sides of 
bulk elements. In this section, it is shown how the integrals over the global domain may be mapped onto 
a one-dimensional master element for evaluation. In particular, we consider an element e whose side, a line 
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element sie, with local nodes 2, 6, 3 forms part of the free surface e G iVi and the extension to other surfaces 
is straightforward. 

Integrals over a line element sie of the liquid-gas free surface are transformed onto a one-dimensional 
master element ^ e [^1, 1] using 



f dsie = / / 



dsie 



^ = ^Jr?+z?, (Ill) 



where the position in the computational domain (ri(^), zi(^)) may be calculated from (77) and derivatives 



with respect to ^, denoted with a prime, are easily found from (103 1. Derivatives of functions with respect 
to sie may then be transformed using 

^ d ' ^ (112) 



dsie dsi, de ^r[^ + z'^ d^ 
The tangent and normal vectors in a free surface element will be 

tl = (tir.^lz) = ^^,2 ^'^2)1/2 ' ni = {nir,ni^) =±{-tiy_,tir), (113) 

where the normal vector is defined to point into the liquid. 

Given this information, the tensor (I — niiii) on the free surface can be determined in terms of the 
master element's coordinates and hence the surface gradient operator can also be obtained: 

(I - nn) = titi + ne.e., = (I - mm) • V ^ ^=^| + (114) 

In a completely analogous way, the expressions required for the evaluation of integrals on the liquid-solid 
interface may be formed. 

8.4.3. Determining local 'line' integrals 

In our computational domain, the contact line corresponds to one point. The vector m, which is defined 
to be normal to the contact line and inwardly tangent to a given surface, will then coincide with the tangent 
vector t of the surface in the (r, 0)-plane. The value of the normal and tangent vectors at the contact line 
may easily be determined by setting the one-dimensional master element's coordinate to the value which it 
will take at the node corresponding to the contact line, i.e. at local node 2 we have ^ = = — 1. 

8.4.4- Evaluating the integrals over the master element 

We have arrived at a set of algebraic equations whose coefficients are integrals over the master element 
which may be easily calculated using numerical integration. The integrands contain products of interpolating 
functions, which are low order polynomials, and hence Gaussian quadrature is well-suited to the task. 
Gaussian quadrature allows the integral to be approximated by a weighted sum of the integrand evaluated 
at a finite number of points, called the Gauss points. Using enough points, one can evaluate polynomial 
expressions of any degree exactly. 

The rules for integration over a one-dimensional master element, say ^ G [—1, 1], are easiest and hence 
presented first. Using Niq Gauss points one can exactly integrate a polynomial of degree 2Niq — 1 using 

(115) 



„1 iViG 



where are the specially chosen Gauss points and Wi are the appropriate weights. Taking Nig — 4 with 

^1 = -0.86113 63115 94052, ^3 = -0.33998 10435 84856, C2 = ^Ci, ^4 -6, 
T^i = 0.34785 48451 37453, W3 = 0.65214 51548 62546, W2 = Wi, Wi = W3, 
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is sufficient for exact integration of all surface integrals. 
For integration over the master element's area we have 



/ / /(r(e,r?),z(e,ry))dCdr?= ^/(6,r/,)W^.. (116) 



One can use tensor products of the one-dimensional case, i.e. integrate over one coordinate and then the 
other, and hence ensure that with Nq Gauss points we can exactly integrate a polynomial of degree 2Nq — 1. 
However, these expressions are not optimal: one can usually integrate a polynomial of degree 2Ng — 1 with 
less than Nq points [73]. The disadvantage of the non-tensor expressions is that the determination of the 
points and the weights becomes more complex. Here, we simply list the Nq = 9 points and weights which 
will exactly integrate polynomials of order 5 and were found to be sufficient for a range of parameter values. 
They are: 





0.00000 00000 00000, 


C4 


0.77459 66692 41483, 






m 


-0.88729 83346 20741, 


V2 


= -0.50000 00000 00000, 


m 


= -0.11270 16653 79258, 


m 


= -0.97459 66692 41483, 


V6 


= -0.80000 00000 00000, 




0.57459 66692 41483, 


Wi 


0.24691 35802 46913, 


W2 


0.39506 17283 95061, 


W4 


0.03478 44646 23227, 




0.05565 51433 97164, 


W7 


0.27385 75106 85414, 




0.43817 20170 96662, 




= Cl: ^5,6 — —£,7.8,9 = £4 


V5 


= m, m = V6, m = va, 


W3 


= Wi, We = Wi, Wq = W'j 



Polynomials of higher order than 5 may exist in the integrands of the bulk integrals, but for a range of 
parameter values the order of quadrature suggested here has been sufficient. It is recommended that a user 
tests different numbers of Gauss points to ensure optimal speed without loss of accuracy. 

When very small elements are used, roughly for element sizes of Imin < 10~®, extreme care has to be 
taken to ensure that the finite precision arithmetic does not create significant errors in the values of the 



integrals. In particular, it was found that for the calculation of the determinant of the Jacobian Jg in (106), 
which is of order Z^i„, and only depends on the relative positions of the nodes in the global domain, the 
coordinate system should be shifted so that the contact line, near which the smallest elements reside, is at 
(0,0). 

8.5. A Newton-Raphson solution method 

At each global node we have a number of algebraic equations whose unknowns are the functions' nodal 
values. These equations and unknowns must each be given a unique position and assembled into a set of 
(nonlinear) algebraic equations to which standard solution methods can be applied. A solution vector a, 
which contains all of the functions' nodal unknowns, can be easily formed by looping through every global 
node and assigning the functions' unknowns there consecutive positions in a, i.e. each unknown in the 
problem is assigned a unique position in a. Then, the residual which determines a given unknown is given 
the same position in a vector containing all of the residual values R. 

To solve the assembled set of (nonlinear) algebraic equations the Newton-Raphson method is used. This 
method takes an approximate solution a™ with associated residuals R{a"^) and then provides an updated 
solution cx™^^ for each entry in the solution vector by solving 

J - c™] = -R{c-), J, = ^1^. (117) 

The method has quadratic convergence and may be repeated until a sufhcient degree of accuracy is achieved. 
In our computations, the resulting linearized equations were solved using the MA41 solver provided by the 
Harwell Subroutine Library. The main complexity of the Newton-Raphson method is in the formation of 



the Jacobian Jy in (117), and this is now described. It may not be necessary to do this at every iteration: 
methods which do not reform the entire Jacobian after each iteration are known as quasi-Newton methods, 
see [801. 
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8.5.1. Formation of the J acobian 

Entries of the Jacobian may be calculated numerically, using a difference formula, or determined ana- 
lytically. Analytic calculation is computationally significantly faster and also improves the convergence of 
the code, particularly when using small elements where the combination of truncation and round-off errors 
in the calculation of the Jacobian can become significant if calculated numerically. Analytic calculation is 
easily achieved for the dependence of the residuals on everything except the free surface unknowns which 
will be considered later. 

Entries of the Jacobian are built, as the residuals are, element-by-element. Below, as a representative 
example, the derivative of the nonlinear local term Ai^Uk = {aikiiui + aiki2Wi)uk, appearing in the r- 
momentum residual (95), with respect to velocity Uj is determined: 



d {Aik{u,w)uk) duk dA,k < c , a /i , i, i « 

^ Aikj, h Uk^ = Aikdjk + UkOikiidji = Aij + UkUikji k = 1-b. (118) 



auj ouj ouj 



As with the residuals, the local nodes will need to be related to the global ones in (117) using the 

function I{e,i). 

The only non-analytic calculation in our code is the dependence of the bulk nodal positions on the free 
surface nodal unknowns, which can become extremely complex when using the bipolar coordinate system 
described in Section [8.2| As the residuals and the Jacobian are evaluated in an element-by-element manner, 
it is most useful to know the dependency of the positions of an element's nodes (r^, z^) {i = 1-6) on the 
global unknowns hj (j = 1-A^i). A quick way of calculating them numerically is by noting the unknowns 
which the nodal positions depend on, in our case the contact line position, apex height and nodal position at 
the end of the spine on which they lie and possibly the adjacent ones. The original nodal positions (r^, Zi) are 
updated by changing the free surface unknown /ij, on which the dependency is required, to hj ^ hj + Ahj, 
where Ahj/hj <^ 1, and obtaining a new mesh with nodal positions (fi,Zj). Repeating this procedure for 
each global nodal unknown (j = l~iVi) gives the required numerically calculated dependencies 

Ti Tj dzi Zi Zi ^119) 



dhj Ahj ' dhj Ah 



and it is noted that (r^, z^) is a local position whilst hj is a global unknown. 

As a representative example, we will show how to calculate the dependence of the pressure terms in the 



z- momentum residual ( 96 ) , Cfj^pk , on a global free surface unknown hj . A possible issue is that the domain 
for the integral depends on the position of the nodes themselves; however, by mapping our integral to the 
fixed master element this problem is resolved. The calculation required is 

where 



Using (110), gives 



dz 



/.fe ^V-fc^r" det Je. (121) 
det Je = ~T,mrm in4>iT J l,m = 1-6 (122) 



so that 



ik 



„, = -^k [Tt„i {ri4>iY + n riTii4>,ri] -ttt-, I, m = 1-6. (123) 
ohj ahj 

It is clear that the derivation and analysis of these integrals is rather cumbersome and prone to error. 
However, one can easily check the correctness of the semi-analytically determined Jacobian by comparing it 
to an approximate one calculated using a forward difference for all entries. This procedure will allow errors 
to be quickly found and rectified. In the FEM, many equations are only dependent on unknowns in adjacent 
elements, so that most entries in the Jacobian are zero. Then, if the equations and unknowns are ordered 
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in a sensible manner, the matrix is sparse and banded with the exception of the dependence of all nodes on 
the contact line position and apex height. 

In some cases, to ensure convergence of the solution, the code was started with the solid moving slowly, 
i.e. |U • e^l ^ 1. Then, using the previous solution as an initial guess, the procedure was repeated, with a 
higher speed, until |U • e^j = 1. This was the case at higher capillary numbers where the final free surface 
shape differs considerably from the initial guess of the static shape. 



.6. Additional details 



From the local asymptotics in ^5.2 it can be seen that both the pressure and normal stress on the solid 
are logarithmically singular as the contact line is approached, and so, strictly speaking, special singular 
elements [HI [S^ should be used to capture this behaviour near the contact line. Such an approach has 
previously been considered for this class of flows in |26j . The singular elements have been incorporated into 
our code, see [31] ; however, to simplify our exposition we have opted against using them here and, as can be 
seen from Figure [jj the agreement of the (integrably singular) asymptotic result for A with the computed 
one is still excellent right up to the node adjacent to the contact line. 
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